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ABSTRACT 

The  lack  of  understanding  of  the  physics  of  the  transition  in  supersonic 
boundary  layer  flows  is  a  major  impediment  in  developing  reliable  transition 
prediction  methods.  This,  in  turn,  represents  a  major  obstacle  in  designing, 
developing,  and  operating  high-speed  flight  vehicles  such  as  the  formerly  pro¬ 
posed  National  Aerospace  Plane  (NASP),  the  future  high-speed  civil  trans¬ 
port  (HSCT),  high-speed  missiles,  high-speed  reconnaissance  aircraft,  and 
the  Theater  Missile  Defense  (TMD)  interceptors.  Transition  to  turbulence 
in  high-speed  flows  is  associated  with  significant  increases  in  aerothermal 
loads  on  the  flight  vehicle,  requiring  additional  thermal  protection.  Super¬ 
sonic  transition  research  suffers  from  the  fact  that  reliable  experiments  are 
difficult  and  very  expensive.  Therefore,  practically  nothing  is  known  about 
the  later,  nonlinear  stages  of  transition  and,  in  particular,  about  the  final 
breakdown  process.  However,  if  and  how  the  final  breakdown  occurs  has 
to  be  understood  in  order  to  identify  which  disturbance  scenarios  lead  to 
breakdown  and  which  do  not.  This  knowledge  is  essential  for  developing  re¬ 
liable  transition  prediction  methods  and  it  is  critical  for  transition  control, 
that  is,  to  delay  transition  and  thus  reduce  the  aerothermal  loads.  To  this 
end,  Direct  Numerical  Simulations  (DNS)  of  supersonic  flow  over  a  flat-plate 
with  and  without  adverse  pressure-gradient  at  Mach  3  were  carried  out  in 
close  collaboration  with  the  experimental  effort  at  Princeton  University  by 
G.  Brown  and  co-workers.  To  confirm  that  simulations  and  experiments 
were  based  on  the  same  “baseflow,”  the  experimental  baseflow  profiles  were 
compared  with  our  Navier— Stokes  results.  The  downstream  development  and 
the  spatial  growth  rates  of  the  disturbances  obtained  from  the  Navier-Stokes 
computations  and  from  experimental  measurements  were  compared  as  well. 
Overall,  a  remarkably  good  agreement  was  achieved.  Towards  the  under¬ 
standing  of  the  nonlinear  mechanisms,  we  investigated  numerous  nonlinear 
resonance  and  breakdown  scenarios.  Our  simulations  have  shown  that  due  to 
the  stabilizing  effects  of  compressibility  for  supersonic  boundary  layers,  the 
transition  process  can  be  stretched  significantly  in  the  downstream  direction 
and  sometimes  the  transition  process  may  even  be  aborted  so  that  a  turbu¬ 
lent  boundary  layer  is  never  fully  established.  The  extent  of  the  transition 
process  and  the  intensity  of  the  temperature  fluctuations,  and  the  resulting 
heat  load,  depend  strongly  on  the  nonlinear  mechanisms. 
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1.  INTRODUCTION 


Transition  in  boundary  layers  at  supersonic  and  hypersonic  speeds  is  a  ma¬ 
jor  unresolved  topic  in  Fluid  Dynamics.  Although  significant  progress  has 
been  made  in  recent  years,  crucial  aspects  of  the  transition  physics  are  still 
in  the  dark.  In  fact,  the  lack  of  a  basic  understanding  of  high-speed  transi¬ 
tion,  and  as  a  consequence  a  lack  of  reliable  transition  prediction  methods, 
was  one  of  the  key  factors  in  the  Defense  Science  Board’s  (1992)  decision 
not  to  recommend  the  construction  of  a  demonstrator  vehicle  for  the  Na¬ 
tional  Aerospace  Plane.  For  the  future  High-Speed  Civil  Transport  (HSCT) 
(Parikh  &  Nagel,  1990),  as  well  as  for  numerous  defense-related  applications 
such  as  high-speed  missiles  (Hingst,  1990;  Korejwo  k  Holden,  1992),  high¬ 
speed  reconnaissance  aircraft  (Scott,  1996),  and  the  Theater  Missile  Defense 
(TMD)  interceptors  (Johnson  et  al. ,  1997),  considerable  progress  toward  the 
understanding  of  high-speed  boundary  layer  transition  is  required  in  order  to 
develop  reliable  transition  prediction  methods  that  can  be  used  for  the  design 
and  safe  operation  of  such  advanced  flight  vehicles.  The  crucial  need  for  re¬ 
liable  transition  prediction  methods  for  high-speed  applications  is  due  to  the 
fact  that  transition  to  turbulence  in  supersonic/hypersonic  boundary  layers 
is  associated  with  considerable  increases  in  heat  transfer.  The  increased  heat 
loads  (caused  by  transition)  on  the  structure  of  the  flight  vehicles  represent 
the  main  difficulties  in  designing  and  operating  high-speed  vehicles.  Appro¬ 
priate  measures  to  guard  against  the  heat  transfer  due  to  aerothermal  loads 
are  expensive  and/or  result  in  significant  weight  penalties.  Good  estimates 
of  the  transition  location  are  of  vital  importance  because  only  then  can  the 
aerothermal  loads  and  surface  temperatures  be  adequately  predicted.  In  ad¬ 
dition  to  surface  heating,  transition  to  turbulence  also  has  a  significant  effect 
on  the  aerodynamic  performance  of  high-speed  flight  vehicles  as  the  skin  fric¬ 
tion  for  turbulent  boundary  layers  is  considerably  higher  than  for  the  laminar 
boundary  layer. 

The  understanding  of  transition  for  low-speed  (incompressible)  boundary 
layers  . is  considerably  ahead  of  that  for  high-speed  (compressible)  boundary 
layers,  although  many  crucial  aspects  are  also  still  not  understood  even  for 
the  low-speed  case.  There  are  several  important  reasons  for  the  considerable 
gap  in  understanding  of  high-speed  transition  relative  to  low-speed  transition. 
Of  course,  historically,  high-speed  flight,  in  particular  hypersonic  flight,  has 
not  been  considered  until  recently  and  therefore  the  need  to  understand  and 
predict  transition  did  not  exist  earlier.  However,  there  are  two  other  main 
reasons  why  it  is  more  difficult  to  obtain  knowledge  for  high-speed  bound¬ 
ary  layer  transition  than  for  the  low-speed  case:  i)  Quality  experiments  for 
high-speed  transition  are  considerably  more  difficult  to  carry  out  than  for 


incompressible  transition  and  require  high-speed  testing  facilities  that  are 
expensive  to  construct  and  expensive  to  operate,  ii)  The  physics  of  high¬ 
speed  boundary  layer  transition  are  much  more  complex  than  for  low  speeds. 
From  linear  stability  theory  (Mack,  1969),  it  is  known  that  multiple  instabil¬ 
ity  modes  exist  for  high-speed  boundary  layer  flows,  in  contrast  to  only  one 
mode  (Tollmien-Schlichting,  TS)  for  the  incompressible  case.  The  so-called 
first  mode  in  supersonic  boundary  layers  is  equivalent  to  the  TS  mode  in  in¬ 
compressible  boundary  layers.  However,  in  contrast  to  incompressible  bound¬ 
ary  layers,  where,  according  to  the  Squire’s  theorem,  two-dimensional  waves 
are  more  amplified  than  three-dimensional  waves,  for  supersonic  bound¬ 
ary  layers  three-dimensional  (oblique)  waves  are  more  amplified  than  two- 
dimensional  ones.  Thus,  experiments  and  theory  always  have  to  deal  with 
the  more  complicated  problem  of  three-dimensional  wave  propagation.  In 
addition  to  the  first  mode,  which  is  viscous,  higher  modes  exist  for  super¬ 
sonic  boundary  layers  that  result  from  an  inviscid  instability  mechanism. 
According  to  linear  stability  theory  (LST),  the  most  unstable  higher  modes 
are  two-dimensional  and  not  oblique,  as  is  the  first  mode.  Also  from  linear 
stability  theory,  it  is  known  that  the  first  mode  is  dominant  (higher  amplifica¬ 
tion  rates)  for  low  supersonic  Mach  numbers  while  for  Mach  numbers  above  4 
the  second  mode  is  dominant  (most  amplified).  In  addition,  for  typical  super¬ 
sonic/hypersonic  flight  vehicle  configurations,  the  three-dimensional  nature 
of  the  boundary  layers  that  occur,  for  example,  on  swept  wings  and/or  lifting 
bodies,  can  give  rise  to  so-called  cross-flow  instabilities  and,  as  a  consequence, 
cross-flow  vortices  that  can  be  stationary  or  traveling. 

Due  to  the  difficulties  in  carrying  out  experiments  (and  “controlled”  experi¬ 
ments,  in  particular)  and  due  to  the  existence  of  multiple  instability  modes, 
the  role  and  importance  of  the  various  instability  modes  in  a  realistic  tran¬ 
sition  process  are  not  understood  at  all.  Of  course,  when  amplitudes  of  the 
various  instability  modes  reach  high  enough  levels,  nonlinear  interactions  of 
these  modes  can  occur.  As  a  consequence,  the  transition  process  in  high¬ 
speed  boundary  layers  is  highly  non-unique  (and  our  simulations  support 
this  conjecture,  see  §  6),  which  means  that  slight  changes  in  the  environment 
or  vehicle  geometry  may  significantly  alter  the  transition  process. 

An  additional  difficulty  arises  from  the  fact  that  for  high-speed  boundary  lay¬ 
ers  the  transition  processes  in  free  flight  may  be  very  different  from  those  in 
the  laboratory.  As  shown  by  Eissler  &  Bestek  (1996),  the  difference  between 
conditions  for  free  flight  (“hot,”  atmospheric  conditions)  and  the  labora¬ 
tory  (“cold,”  laboratory  conditions)  has  a  considerable  effect  on  the  stability 
behavior  and,  as  a  consequence,  on  the  transition  processes.  This  is  best 
summarized  by  a  quote  from  Stetson  (1990),  a  pioneer  in  experimental  high¬ 
speed  transition  research:  . .  one  should  not  expect  a  transition  Reynolds 
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number  obtained  in  any  wind  tunnel,  conventional  or  quiet,  to  be  directly 
relatable  to  flight.”  This  fact  clearly  defines  the  usefulness  and  critical  im¬ 
portance  of  numerical  simulations.  The  numerical  simulation  codes  can  be 
tested  and  validated  by  detailed  comparison  with  laboratory  experiments. 
Thereafter,  they  can  then  be  applied  with  confidence  to  predict  the  effects 
of  free- flight  conditions  on  the  transition  processes  and  the  resulting  aero¬ 
dynamic  and  aero-  thermodynamic  behaviors.  Thus,  only  simulations  can 
provide  the  crucial  understanding  and  information  for  design  and  safe  oper¬ 
ation  of  high-speed  vehicles. 


2.  RELATED  PREVIOUS  RESEARCH 


2.1  Theoretical  Investigations 

There  are  a  large  number  of  scientific  publications  available  on  transition 
research,  with  the  largest  number  of  them  dealing  with  low-speed  transition. 
Recent  investigations  of  transition,  both  high-  and  low-speed,  were  presented 
at  the  IUTAM  Symposium  on  Laminar  Turbulent  Transition  (Fasel  &  Saric, 
2000).  Some  of  the  most  important  aspects  of  high-speed  transition  that 
are  relevant  to  the  proposed  research  are  discussed  below.  The  main  body 
of  the  present  knowledge  concerning  high-speed  transition  is  still  what  was 
obtained  using  Linear  Stability  Theory  (LST)  by  Mack  (1969,  1975,  1984, 
2000).  Based  on  the  findings  by  Mack,  the  linear  stability  behavior  of  com¬ 
pressible  (supersonic)  boundary  layers  differs  from  the  incompressible  case 
in  several  significant  aspects: 

i.  More  than  one  instability  mode  exists  for  Ma  >  2.2:  the  first  mode  and 
the  second  and  higher  (multiple)  modes. 

ii.  The  first- mode  disturbances  are  viscous  (vortical)  and  are  similar  to  the 
Tollmien-Schlichting  (TS)  modes  of  incompressible  boundary  layers. 
First-mode  disturbances  dominate  (largest  amplification  rates)  at  low 
supersonic  Mach  numbers.  However,  in  contrast  to  the  incompressible 
case,  the  most  amplified  first-mode  disturbances  are  three-dimensional 
(oblique)  and  not  two-dimensional. 

iii.  The  second  and  higher  modes  are  inviscid  (acoustic)  and  dominate  at 
Mach  numbers  higher  than  about  4.  The  most  unstable  modes  here 
are  always  two-dimensional  (in  contrast  to  the  first  mode). 

iv.  In  addition  to  the  inviscid  (acoustic)  higher  modes,  Mack  identified 
additional  viscous  modes  (Viscous  multiple  solutions’)  which,  to  date, 
have  not  been  identified  in  experiments.  However,  they  were  also  found 
in  the  Direct  Numerical  Simulations  (DNS)  of  Eissler  &  Bestek  (1996). 

v.  First-mode  disturbances  can  be  attenuated  (as  for  the  incompressible 
case  for  air)  by  wall  cooling,  wall  suction,  and  favorable  pressure  gra¬ 
dients  (Malik,  1989). 

vi.  The  second  and  higher  inviscid  modes  can  be  stabilized  by  favorable 
pressure  gradients  and  suction;  however,  they  are  destabilized  by  wall 
cooling. 


For  linear  stability  theory  analysis,  the  effects  of  the  growing  boundary  layer 
on  the  disturbance  growth  are  typically  neglected  (“parallel  theory”).  How¬ 
ever,  non-parallel  effects  can  be  included  by  using  the  Parabolized  Stability 
Equations  (PSE)  approach  (Bertolotti,  1991;  Bertolotti  et  al.,  1992;  Herbert, 
1994).  Depending  on  various  parameters  (Mach  number,  Reynolds  number, 
frequency,  etc.),  non-parallel  effects  can  significantly  influence  the  distur¬ 
bance  growth  rates.  LST  and  linear  PSE  are  only  applicable  for  the  first 
(linear)  stage  of  the  transition  process  where  disturbance  amplitudes  are 
small  and  nonlinear  interactions  are  negligible.  Nonlinear  PSE,  on  the  other 
hand,  is  applicable  to  the  nonlinear  stages  of  transition  (Bertolotti  et  al., 
1992;  Herbert,  1994),  although  the  computational  effort  becomes  overwhelm¬ 
ing  when  the  development  becomes  strongly  nonlinear  in  the  later  stages  of 
transition.  Also,  analogous  to  incompressible  boundary  layer  transition,  sev¬ 
eral  attempts  have  been  made  to  apply  secondary  instability  theory  to  model 
the  initial  three-dimensional  nonlinear  development  (see,  Masad  k  Nayfeh, 
1990;  El-Hady,  1991,  1992;  Ng  k  Erlebacher,  1992,  for  example).  However, 
whether  or  not  any  of  these  secondary  instability  mechanisms  are  relevant 
for  supersonic  transition  is  an  open  question  because  it  is  very  difficult  to 
unequivocally  identify  them  in  experiments  (see  §  2.2). 

2.2  Experimental  Investigations 

Transition  experiments  in  high-speed  flows  are  extremely  difficult  and  very 
expensive.  Therefore,  relatively  few  successful  experimental  efforts  have  been 
discussed  in  the  open  literature.  Most  experiments  have  emphasized  inves¬ 
tigations  of  the  linear  regime  and  the  early  stages  of  the  transition  process. 
Some  examples  are  the  experiments  by  Laufer  k  Vrebalovich  (1960);  Kendall 
(1975);  Stetson  et  al.  (1983);  Stetson  (1988);  Kosinov  et  al.  (1990);  Stetson 
k  Kimmel  (1992a);  Schneider  et  al.  (1996),  and  Corke  (see  Cavalieri,  1995). 
An  overview  of  the  experimental  efforts  up  to  1992  is  given  by  Stetson  k 
Kimmel  (19926).  The  experiments  essentially  verified  some  of  the  important 
parts  of  the  linear  theory.  However,  quantitative  differences  often  occur  that 
may  be  explained  by  the  fact  that  in  the  experiments  the  transition  pro¬ 
cess  was  “natural,”  i.e.,  it  was  initiated  by  the  environmental  disturbances, 
and  not  by  “controlled”  disturbance  input  (analogous  to  a  vibrating  ribbon 
as  in  incompressible  transition  experiments).  Also,  quantitative  differences 
between  experimental  measurements  and  LST  may  be  caused  by  the  effects 
of  the  growing  boundary  layer  being  neglected  in  the  linear  stability  theory 
analysis  (“parallel  theory”). 

All  experimental  efforts  suffered  more  or  less  from  difficulties  in  control¬ 
ling  the  disturbance  environment  (such  as  sound  radiated  from  turbulent 


boundary  layers  on  the  tunnel  walls).  Nevertheless,  these  “natural”  transi¬ 
tion  experiments  could  identify  first  and  second  instability  modes  (Kendall, 
1975;  Stetson  et  al. ,  1983).  However,  for  example  for  Mach  2,  considerable 
discrepancies  arose  between  boundary  layers  on  a  flat  plate  and  boundary 
layers  on  axisymmetric  cones  when  “blow  down”  facilities  were  used.  For  ax¬ 
isymmetric  cones,  high-frequency  second  modes  were  dominant,  while  only 
low-frequency  first  mode  disturbances  were  observed  for  the  flat  plate  (Stet¬ 
son  &  Kimmel,  19926).  In  contrast,  in  an  experiment  using  a  Ludwieg  tube 
for  a  sharp-nosed  cone  at  Mach  5,  no  dominant  second-mode  disturbances 
could  be  detected  (Wendt,  1993). 

With  the  more  recent  experiments  for  a  flat  plate  and  axisymmetric  cones 
at  Mach  number  3.5  in  the  NASA-Langley  “Quiet  Tunnel,”  a  number  of  dis¬ 
crepancies  between  LST  and  other  experiments  were  resolved  (Chen  et  al, 
1989;  Cavalieri,  1995).  Indications  of  nonlinear  developments  in  the  transi¬ 
tion  process  were  observed  by  Stetson  et  al.  (1983)  for  a  cone  at  Mach  8. 
Most  of  the  experimental  efforts  suffered  from  the  deficiency  that  no  “con¬ 
trolled”  disturbances  could  be  introduced  to  allow  for  detailed  quantitative 
comparisons  with  linear  theory  and,  in  particular,  to  allow  for  systematic 
investigations  of  the  nonlinear  stages  of  transition.  Because  of  the  lack  of  ex¬ 
perimental  evidence  concerning  the  process  in  the  later  stages  of  transition, 
it  is  totally  unclear  which  instability  modes  and  which  nonlinear  mechanisms 
are  responsible  for  the  final  breakdown  to  turbulence  in  supersonic  boundary 
layers.  The  controlled  experiments  for  a  Mach  2  boundary  layer  by  Kosinov 
&  Tumin  (1996)  using  a  harmonic  point  source  for  the  disturbance  excitation 
have  indicated,  however,  that  secondary  instability  mechanisms  were  present. 
In  fact,  Kosinov  &  Tumin  speculated  that  it  was  a  subharmonic  resonance 
with  “oblique”  fundamental  disturbances.  More  recently,  Maslov  and  co¬ 
workers  (see  Shiplyuk  et  al.,  2003)  reported  “controlled”  experiments  for  a 
sharp-nosed  cone  at  Ma  =  5.95  using  a  glow-  discharge  actuator  to  generate 
harmonic  point  source  disturbances.  They  investigated  several  nonlinear  in¬ 
teractions  and  identified  a  “classical”  subharmonic  resonance  (with  the  2-D 
second  mode  as  primary  disturbance)  as  a  possible  breakdown  mechanism 
involving  possibly  a  3-D  first  mode  as  the  subharmonic.  However,  in  order 
to  confirm  this  conjecture,  detailed  spatial  and  temporal  resolution  of  the 
measurements  would  be  required  which,  of  course,  is  difficult  experimentally. 
Shiplyuk  et  al.  (2003)  state  in  their  paper,  that  numerical  calculations  would 
be  helpful  to  clarify  the  scenarios  of  nonlinear  interactions  that  are  identified 
in  the  present  work. 


2.3  Numerical  Simulations 


Due  to  the  difficulties  in  experimental  investigations  of  high-speed  boundary- 
layer  transition  and  due  to  the  limitations  of  linear  stability  theory,  so-called 
Direct  Numerical  Simulations  (DNS)  represent  a  promising  tool  for  high¬ 
speed  transition  research.  In  DNS  the  complete  Navier-Stokes  equations  are 
solved  directly  by  proper  numerical  methods  without  making  restricting  as¬ 
sumptions  with  regard  to  the  baseflow  and  the  form  and  amplitude  of  the 
disturbance  waves.  Therefore,  DNS  is  particularly  well  suited  for  investiga¬ 
tions  of  the  nonlinear  development  that  is  characteristic  of  the  later  stages  of 
high-speed  boundary  layer  transition.  Two  fundamentally  different  models 
are  used  for  DNS:  First,  the  so-called  “temporal  model,”  which  is  based  on  the 
assumption  that  the  baseflow  does  not  change  in  the  downstream  direction 
(thus  excluding  nonparallel  effects).  Also,  assuming  spatial  (downstream) 
periodicity  of  the  disturbances,  the  disturbance  development  (growth  or  de¬ 
cay)  is  then  in  the  time-direction.  The  temporal  model  is  analogous  to  the 
temporal  approach  in  LST  with  the  frequency  being  complex  and  the  spatial 
wave  number  real.  Due  to  the  underlying  assumptions,  the  temporal  model 
can  only  provide  qualitative  results.  On  the  other  hand,  since  a  relatively 
short  integration  domain  can  be  used  in  the  downstream  direction  (typically 
one  or  two  wave  lengths  of  the  fundamental  wave)  temporal  simulations  are 
relatively  inexpensive. 

In  contrast,  in  the  “spatial  model”  no  assumptions  are  made  with  regard 
to  the  baseflow  (thus  nonparallel  effects  are  included).  The  disturbance  de¬ 
velopment  (growth  or  decay)  is  in  the  downstream  direction  as  in  physical 
laboratory  or  free-flight  situations.  Thus,  the  spatial  model  allows  realistic 
simulations  of  high-speed  transition  and  direct  comparison  with  wind-tunnel 
or  free-flight  experiments.  However,  simulations  based  on  the  spatial  model 
are  typically  much  more  costly  than  for  the  temporal  model  because  a  much 
larger  downstream  integration  domain  is  required  (many  wave  lengths  of  the 
fundamental  disturbance  wave).  This  is  particularly  true  for  simulations  of 
high-speed  boundary  layer  transition,  where  the  growth  rates  of  the  distur¬ 
bance  waves  are  much  smaller  than  for  the  incompressible  case  and  where 
the  growth  rates  of  certain  modes  decrease  with  increasing  Mach  number. 
Thus,  relatively  large  (in  the  downstream  direction)  integration  domains  are 
required  to  allow  small  disturbances  to  grow  to  the  large  amplitudes  that 
characterize  the  nonlinear  stages  of  the  transition  process  and  finally  lead  to 
the  breakdown  to  turbulence.  As  a  consequence,  spatial  simulations  of  high¬ 
speed  transition  are  computationally  very  challenging.  Detailed  discussions 
of  the  DNS  methodology  for  investigations  of  boundary  layer  transition,  in 
particular,  discussions  of  the  temporal  and  spatial  approach,  are  given  by 


Fasel  k  Konzelmann  (1990)  ,  Kleiser  k  Zang  (1991),  and  Reed  (1993). 
Probably  the  first  transition  simulation  for  supersonic  boundary  layers,  al¬ 
though  restricted  to  two-  dimensional,  yet  spatially  evolving,  disturbances, 
was  by  Bayliss  et  al.  (1985),  who  employed  an  approach  analogous  to  that  by 
Fasel  (1976)  for  incompressible  boundary  layers.  The  first  three-dimensional 
yet  temporal  DNS  for  flat-plate  high-speed  boundary  layer  transition  was 
performed  by  Erlebacher  k  Hussaini  (1990).  Here,  only  the  linear  and  early 
nonlinear  stages  were  explored.  Other  temporal  simulations  were  performed 
by  Normand  &  Lesieur  (1992),  Pruett  &  Zang  (1992),  Dinavahi  k  Pruett 
(1993),  Adams  &  Kleiser  (1993).  From  such  temporal  simulations,  Normand 
k  Lesieur  (1992)  found  that,  for  their  case  of  a  fiat-plate  boundary  layer  with 
Ma  =  5,  transition  occurred  via  a  subharmonic  secondary  instability  for  the 
second  mode.  This  finding  was  consistent  with  results  from  simulations  by 
Adams  k  Kleiser  (1993),  Pruett  &  Zang  (1992),  and  Dinavahi  k  Pruett 
(1993)  for  a  boundary  layer  at  Mach  4.5  on  a  hollow  cylinder  (the  axisym- 
metric  analog  of  a  flat- plate  boundary  layer).  However,  the  main  weakness  of 
these  “temporal”  simulations  is  the  fact  that  they  do  not  take  the  boundary 
layer  growth  into  account.  In  fact,  experiments  by  Stetson  k,  Kimmel  (1993) 
and  PSE  calculations  by  Chang  et  al.  (1991)  indicated  that  subharmonic 
resonance  may  not  be  the  preferred  route  to  transition  in  realistic,  growing 
boundary  layers  (which  include  nonparallel  effects). 

Realistic  simulations  of  transition  scenarios  including  the  effects  of  the  grow¬ 
ing  boundary  layer  require  the  use  of  the  spatial  simulation  model.  The 
first  three-dimensional  spatial  simulations  of  transition  in  supersonic  bound¬ 
ary  layers  were  reported  by  Thumm  (1991)  for  a  Mach  number  of  1.6.  In 
fact,  from  these,  and  follow-up  simulations  (Fasel  et  al.,  1993),  it  was  dis¬ 
covered  that  a  new  “Oblique  Breakdown”  mechanism  produces  much  larger 
growth  rates  than  either  subharmonic  or  fundamental  resonance  and  requires 
much  lower  disturbance  amplitudes.  Therefore,  the  “Oblique  Breakdown”  is 
a  likely  candidate  for  a  viable  path  to  transition  for  supersonic  boundary  lay¬ 
ers.  Using  PSE  calculations,  Chang  k  Malik  (1993)  confirmed  the  validity  of 
this  oblique  breakdown  for  a  flat-plate  boundary  layer  at  Ma  =  1.6.  Based 
on  our  DNS  code  (see  Fasel  et  al.,  1993),  Bestek  &  Eissler  (1996)  performed 
simulations  for  Mach  4.8  and  investigated  various  nonlinear  mechanisms  in¬ 
cluding  the  “oblique  breakdown”  mechanism.  Bestek  k  Eissler  also  found, 
for  the  first  time,  an  additional  “higher  viscous”  mode,  which  Mack  (1969) 
had  predicted  using  LST  analysis.  Pruett  k  Chang  (1993)  carried  out  spa¬ 
tial  DNS  for  a  flat-plate  boundary  layer  at  Mach  4.5  and  provided  a  detailed 
comparison  with  PSE  results.  Later,  an  improved  version  of  the  code  (Pruett 
et  al.,  1995)  was  applied  to  a  simulation  of  transition  on  axisymmetric  sharp 
cones  at  Mach  8  (Mach  6  after  the  shock;  Pruett  k  Chang,  1995).  The 


simulation  was  combined  with  PSE  calculations  such  that  the  linear  and 
moderately  nonlinear  stages  were  computed  by  PSE  while  the  strongly  non¬ 
linear  and  breakdown  stages  of  transition  were  computed  by  spatial  DNS. 
This  approach  was  motivated  by  the  experience  that  linear  and  moderately 
nonlinear  wave  propagations  can  be  computed  more  efficiently  with  PSE 
while  the  strongly  nonlinear  and  breakdown  stages  (requiring  many  spariwise 
Fourier  modes)  are  computed  more  efficiently  with  DNS.  In  this  simulation,  a 
second-mode-breakdown  resonance  was  also  investigated.  The  so-called  rope¬ 
like  structures  obtained  from  numerical  flow  visualizations  of  the  simulation 
data  for  this  breakdown  process  are  similar  to  those  observed  in  high-speed 
transitional  boundary  layers  on  cones  (see  Pruett  &  Chang,  1995).  Using 
numerical  simulations,  the  leading  edge  receptivity  of  high-speed  boundary 
layers  has  been  investigated  extensively  by  Zhong  and  co-workers  (see  Zhong, 
2001).  More  recently,  the  magnetic  field  effects  on  the  second- mode  insta¬ 
bilities  for  a  weakly  ionized  boundary  layer  at  Ma  =  4.5  using  numerical 
simulations  were  investigated  by  Cheng  et  al.  (2003). 

Relatively  few  attempts  have  been  made  to  employ  Large-Eddy  Simulation 
(LES)  for  transitional  flows  in  supersonic  boundary  layers  (see,  for  example, 
Normand  &;  Lesieur,  1992;  Zang  et  al.,  1992).  These  simulations  reported 
in  the  literature  have  to  be  viewed  as  being  of  an  exploratory  nature  with 
regard  to  the  applicability  of  LES  for  boundary  layer  transition.  Nevertheless, 
these  attempts  demonstrated  that  LES  could  be  employed  advantageously 
for  supersonic  transition  simulations.  The  main  issues  in  applying  LES  for 
supersonic  transition  simulations  are  the  use  of  proper  subgrid-scale  (SGS) 
models  that  are  physically  consistent  through  the  entire  transition  process. 
Of  course,  LES  would  be  a  highly  valuable  tool  for  investigating  the  late 
stages  of  transition.  As  mentioned  previously,  DNS  for  supersonic  boundary 
layer  transition  is  very  expensive  due  to  the  large  computational  domains  that 
are  required  in  the  downstream  direction.  In  addition,  for  simulations  that 
include  the  final  stages  of  transition  (the  actual  breakdown  to  turbulence),  an 
extremely  fine  grid  would  be  required  which,  as  a  consequence,  would  place 
high  demands  with  regard  to  computer  memory  and  computation  times. 
LES  on  the  other  hand  would  require  considerably  less  resolution  and,  as  a 
consequence,  the  amount  of  computer  memory  and  computation  times  would 
be  reduced  accordingly. 


3.  COMPUTATIONAL  METHOD 


The  time-dependent  flow  over  a  flat  plate  is  investigated  using  the  spatial 
DNS  model.  This  model  allows  for  a  direct  comparison  with  experiments 
where  disturbances  develop  in  the  downstream  direction.  The  complete 
Navier-Stokes  equations  for  compressible  flow  form  the  governing  equations 
and  thus  include  all  non-parallel  and  nonlinear  effects.  In  this  chapter,  details 
on  the  governing  equations,  boundary  and  initial  conditions,  buffer  domains 
and  introduction  of  controlled  disturbances  are  given  followed  by  a  brief  sum¬ 
mary  of  the  numerical  method.  Details  on  the  numerical  method  and  code 
performance  can  be  found  in  Harris  (1997)  and  von  Terzi  (2004). 


3.1  Governing  Equations 


Conservation  of  mass  (continuity  equation),  conservation  of  momentum,  and 
conservation  of  total  energy  form  the  set  of  governing  equations.  These 
are  applied  to  a  three-dimensional  Cartesian  coordinate  system.  The  non- 
dimensional  form  of  the  governing  equations  is  obtained  by  by  using  a  length 
scale  (the  length  of  the  plate  in  the  experiments  L*^)  and  free-stream  values 
of  velocity,  temperature,  density,  and  specific  heat  ([/£,,  T^,  p*x  and  C*^, 
respectively)  for  normalization  of  all  variables.  The  non-dimensional  length 
of  the  flat  plate  becomes  equal  to  one  and  the  following  non-dimensional 
variables  are  obtained: 
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The  Mach  and  Reynolds  numbers  are  introduced  as 
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Ma  =  — 22 
a* 


U* 


Re  = 
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Pi 


(3.1) 


respectively.  The  non-dimensional  continuity,  momentum  and  energy  equa¬ 
tions  are  given  by 

dp  d 

-  +  —  ^  =  0 


(3.2) 
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dEt  d 

~dt+dxj  ^Et  +  P^Uj  ~  UiTij  +  =  0  ‘ 

The  total  energy,  viscous  stress  and  heat  flux  terms  are 
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respectively.  The  viscosity  is  obtained  using  Sutherland’s  law  and  the  equa¬ 
tion  of  state  becomes 

_  PT 


3.2  Computational  Domain 


Figure  3.1  Computational  domain  and  boundaries. 


The  equations  are  discretized  and  solved  in  the  computational  domain  shown 
in  figure  3.1.  The  dashed  lines  indicate  the  boundary  layer  that  forms  when 
a  uniform  flow  with  velocity  encounters  the  flat  plate.  An  example  for  a 
computational  domain  is  represented  by  a  box  which  starts  at  Xq  and  ends 
at  .  The  value  for  the  domain  height  is  Dm  ■  Disturbances  are  introduced 
between  the  locations  X\  and  Xo .  All  simulations  performed  in  this  work  deal 
with  3D  waves  generated  in  pairs  in  the  disturbance  slot  with  T  and  — ^  as 


their  wave  angle.  The  wave  angle  determines  their  spanwise  wave  number 
since 

tan(tf)  =  2-  .  (3.9) 

The  wave  with  T  and  the  wave  with  — $  as  its  wave  angle  have  the  same 
spanwise  wave  number  7.  The  domain  width  is  chosen  such  that  it  matches 
half' of  one  wavelength  A*  of  the  3D  waves.  Locations  X3,  x4  and  x$  specify 
the  start  and  end  locations  of.  the  pressure  gradient  (if  applied) .  They  are 
explained  in  more  detail  in  §  3.4.4. 


3.3  Initial  Conditions 

All  simulations  are  carried  out  in  two  parts.  First,  a  “baseflow”  is  computed 
as  the  steady  solution  to  the  Navier-Stokes  equations  and,  second,  this  base- 
flow  is  used  as  an  initial  condition  for  the  stability  investigations,  where 
controlled  disturbances  are  then  introduced.  For  the  (unsteady)  forced  flow 
simulations,  the  baseflow  is  also  required  as  input  for  the  inflow,  outflow  and 
free-stream  boundary  conditions  (see  §  3.4). 


3.3.1  Using  Compressible  Blasius  Similarity  Solution 

All  simulations  at  supersonic  speeds  use  the  compressible  boundary  layer  sim¬ 
ilarity  solution  as  initial  condition  for  obtaining  the  baseflow.  The  similarity 
solution  is  obtained  from  the  Blasius  boundary  layer  ordinary  differential 
equation  (ODE)  where  the  prime  denotes  a  derivative  with  respect  to  the 
independent  (similarity)  variable  77. 

(Cf")'  +  ff"  =  0  (3.10) 


(kS')  +fs'  +  C{ff  =  0  (3.11) 

with 

c  =  £t 1 

Pl oPlo 

/  =  f(v) 

9  =  9{v) 


Several  simplifications  in  deriving  this  ODE  were  made  so  that  the  Navier- 
Stokes  solution  will  differ  from  this  initial  condition.  These  perturbations  will 
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travel  downstream  and,  finally,  will  be  convected  out  of  the  computational 
domain.  The  resulting  flow  field  is  the  desired  baseflow  for  the  simulations 
without  a  pressure-gradient.  However,  for  the  simulations  with  an  adverse 
pressure-gradient  (APG)  in  the  downstream  direction,  the  above  “baseflow” 
is  taken  as  an  initial  guess.  Then,  a  constant  pressure  gradient  was  introduced 
through  the  free-stream  boundary  condition  and  the  baseflow  with  APG  is 
computed  (see  figure  3.2). 


Figure  3.2  Generation  of  a  compressible  baseflow  with  APG. 


3.3.2  Using  Falkner-Skan  Similarity  Solution 

For  the  incompressible  validation  simulation  with  an  adverse  pressure  gra¬ 
dient  in  the  streamwise  direction,  the  solution  to  the  third-order  ODE  of 
Falkner  &  Skan  (1931)  (c.f.  Schlichting  fe  Gersten,  2000)  is  used  as  the  ini¬ 
tial  condition  for  obtaining  the  baseflow.  These  equations  were  derived  by 
introducing  the  similarity  variable 
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as  the  independent  variable  and  f  as  the  dependent  variable  with 
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(again  a  prime  denotes  a  derivative  with  respect  to  rj).  This  transformation 
is  only  possible  if  the  following  assumption  for  the  free-stream  velocity  Ue(x) 
is  made: 

Ue(x)  =  U0xm,  (3.15) 

with  Uo  =  x0m  and  m  =  /?#/ (2  —  /?#)  where  /?#  is  the  Hartree  parameter. 
For  pH  =  —0.1988  the  flow  separates. 


The  rc-momentum  equation  provides  an  expression  for  the  pressure  gradient 


dp 

dx 


(3.16) 


Integrating  this  expression  in  x  leads  to  the  pressure  distribution  at  the  free- 
stream 

with  l/(«Ma2)  as  the  integration  constant.  Equation  (3.17)  is  the  final  equa¬ 
tion  which  is  implemented  into  the  compressible  DNS  code.  The  integration 
constant  represents  the  value  of  the  pressure  at  the  free-stream  at  location 
x  =  x0  and  is  obtained  by  equation  (3.8)  for  p  =  T  =  1.  The  Mach  number 
appears  in  equation  (3.17)  since  the  DNS  code  used  for  the  present  research 
is  compressible. 

Finally,  the  third-order  ODE  derived  by  Falkner  &  Skan  has  the  form 


f'  +  ff  +  M  l-/*)  =  0 


with  the  boundary  conditions 


(3.18) 


/  (0)  =  0;  /'( 0)  =  0;  /'( oo)  =  1.  (3.19) 


It  is  solved  numerically  and  the  solution  serves  as  an  initial  condition  for  the 
DNS  code. 

Since  the  DNS  code  solves  the  compressible  Navier-Stokes  equations  and 
the  initial  condition  is  incompressible,  a  value  for  the  temperature  has  to  be 
computed.  The  total  enthalpy  ht  is  assumed  to  be  constant  over  y.  This 
means  that  the  enthalpy  in  the  free-stream  is  equal  to  the  enthalpy  in  the 
boundary  layer 

h*te  =  h*t  =  const.  (3.20) 

The  substitution  of  h*t  =  c*pT*  +  1/2  (u*2  +  v*2)  into  equation  (3.20)  gives 


c*pT* 


+ 


(u*2  +  V**) 


=  c*T* 


Ue 

-i - - 

2 


*2 


which  is  put  in  non-dimensional  form  and  solved  for  T 


(3.21) 


T  =  1  +  ^-~Ma2(U?  -  u2  -  v s);  (3.22) 

where  v  has  the  form 


v  = 


m  +  1  Ue(x) 
2  x  Re 


fl—m 
\m  +  1 


(3.23) 


Table  3.1  summarizes  the  equations  that  are  used  in  the  DNS  code  to  generate 
an  appropriate  initial  condition  for  the  incompressible  validation  case  with 
adverse  pressure  gradient. 


quantity 

equation 

UJC 

f'Ue 

vie 

/“‘Stei/'-/) 

Pic 

i(l -(£)*”) +  1/(«M>S) 

Tic 

1  +  s=i Ma\U!  -  u)c  -  vja) 

Pic 

K  Ma:  picjTjc 

Table  3.1  Equations  used  for  initial  condition  of  incompressible  validation 
case  with  APG. 


3.4  Boundary  conditions 

The  computational  domain  and  its  boundaries  are  shown  in  figure  3.3.  In 
the  lateral  direction  z,  periodicity  is  assumed.  More  details  about  the  lateral 
boundary  can  be  found  in  von  Terzi  (2004).  The  baseflow  and  the  forced 
flow  simulations  are  computed  with  the  same  boundary  conditions  with  the 
exception  of  the  free-stream,  where  for  the  baseflow  the  pressure  is  fixed  and 
for  the  unsteady  simulations  a  decay  condition  is  employed  for  all  disturbance 
quantities. 


free-stream  boundary 


3.4.1  Inflow 

Incompressible  Flow 

At  the  inflow,  all  quantities  except  for  the  pressure  are  fixed  with  the  value 
of  the  initial  condition 


<p  =  4>Ic. 


(3.24) 


For  the  pressure,  a  gradient  condition  is  employed.  Because  of  the  pressure 
gradient  in  streamwise  direction,  its  derivative  is  not  equal  to  zero.  The 
change  of  the  pressure  in  y  is  neglected  so  that  its  derivative  in  x  has  the 
same  value  as  at  the  free-stream  over  the  entire  inflow 

||  =  -U^mx2™-1  and  ^  =  0 ,  (3.25) 

for  flow  with  and  without  APG,  respectively.  Since  T  and  p  are  fixed  the 
equation  of  state  is  violated.  Nevertheless,  upstream- travehng  sound  waves 
can  pass  through  the  boundary. 


Compressible  Flow 


For  a  compressible  flow  all  quantities  are  fixed  to  the  value  of  the  initial 
condition 


<P  =  <f>ic, 


(3.26) 


except  for  the  subsonic  region  of  the  boundary  layer  where  the  pressure 
derivative  is  set  to  zero 

|=0.  (3.27) 

This  is  not  valid  for  a  streamwise  pressure  gradient  but  the  subsonic  region 
is  very  small  and  the  ill-posed  problem  did  not  produce  any  difficulties. 


3.4.2  Outflow 

Incompressible  Flow 

The  outflow  is  treated  by  applying 

^  =  0.  (3.28) 

where  <p  symbolizes  the  conservative  variables  p,  pUi  and  E.  Since  the  pres¬ 
sure  is  not  held  constant  at  the  inflow,  a  Dirichlet  condition  must  be  employed 
to  the  outflow 

P  =  Pic ,  (3.29) 

otherwise  the  pressure  field  drifts  from  its  specified  value.  For  the  incompress¬ 
ible  validation  case,  the  fixed  pressure  was  responsible  for  strong  reflections 
of  pressure  waves.  The  reflections  caused  upstream-traveling  waves  in  all 
other  quantities  that  were  reflected  at  the  inflow.  A  remedy  for  this  problem 
was  found  by  ramping  down  the  disturbances  before  they  reach  the  inflow  or 
outflow  boundary  (see  §  3.5). 


Compressible  Flow 

For  all  conservative  variables  as  well  as  for  the  pressure,  the  condition 


0V 

dx 2 


(3.30) 


is  applied.  Reflections  of  disturbances  were  very  weak,  however  this  condition 
fails  where  u  <  0  or  u  &  0.  For  that  reason,  a  buffer  domain  (see  §  3.5), 
which  ramps  the  disturbances  to  £ero,  had  to  be  used,  in  particular  for  the 
oblique  breakdown  investigations. 


3.4.3  Wall 

The  incompressible  and  compressible  calculations  are  performed  with  the 
same  wall  boundary  conditions.  The  previously  described  computational 
boundary  conditions  are  artificial  and  usually  do  not  represent  the  physical 
solution  at  the  enforced  location.  An  exception  for  that  is  the  wall.  The 
no-slip  and  no-penetration  conditions  are  used,  thus 


u  =  v  =  w  =  0. 


(3.31) 


Note  that  v  is  non-zero  at  the  disturbance  slot.  The  wall  is  assumed  to  be 
adiabatic,  i.e., 

&T  n 

an  = 1 °'  <3-32) 

and  the  pressure  is  calculated  from  the  y-momentum  equation 


dp 

dy 


d  2  d 
dy(pV  ^  +  dx^Txv)  + 


dy  (TW)  +  dZ^' 


(3.33) 


Note  that  the  time  derivative  is  zero  for  the  steady  baseflow  and  neglected 
at  the  forcing  slot  for  the  simulations  with  forcing. 


3.4.4  Free-Stream 

One  focus  of  this  report  is  the  investigation  of  the  influence  of  an  APG  on 
transition  for  a  compressible  boundary  layer.  For  this  purpose,  an  APG 
has  to  be  introduced.  The  simplest  way  of  achieving  this  goal  is  to  provide 
an  appropriate  free-stream  boundary  condition.  In  experiments  a  pressure 
gradient  can  be  created  at  the  free-stream  through  geometrical  changes.  For 
supersonic  flows,  a  decreasing  of  the  cross-sectional  area  of  the  flow  passage 
and,  for  subsonic  flows,  an  increasing  of  the  cross-sectional  area  of  the  flow 
passage  generate  an  APG. 


Incompressible  Flow 


For  the  baseflow  calculation,  as  an  initial  condition,  a  Falkner-Skan  similarity 
solution  is  employed  (see  §  3.3.2)  which  already  has  the  correct  pressure 
distribution  at  the  free-stream.  Therefore  p  only  needs  to  be  fixed  to  the 
initial  value 

P  =  Pic  •  (3.34) 

For  u  and  p,  von  Neumann  conditions  are  used,  hence 


(3.35) 


and 


(3.36) 


This  is  a  good  approximation  because  the  density  should  be  constant  (in¬ 
compressible  flow)  and  u  should  not  change  in  y  outside  of  the  boundary 
layer.  From  the  continuity  equation  the  derivative  of  the  v- velocity  can  be 
obtained, 

dv  _  du  dUe{x) 

dy  dx  dx 

Using 

Ue(x)=U0xm 

equation  (3.37)  yields 

dv  , 

■x-  =  -U0mxm  1 . 
dy 


(3.37) 


(3.38) 


(3.39) 


For  the  unsteady  simulations,  an  important  issue  is  the  imposition  of  the 
pressure  gradient.  A  boundary  condition  at  the  free-stream  should  not  mod¬ 
ify  the  specified  pressure  distribution  of  the  baseflow  and,  at  the  same  time, 
allow  disturbances  to  pass  the  domain  boundary.  Kloker  (1993),  in  his  code, 
divided  the  flow  quantities  into  a  disturbance  flow  and  a  baseflow.  For  the 
disturbance  quantities  he  applied  a  decay  condition.  The  u- velocity  of  the 
baseflow  was  set  to  the  similarity  value  of  equation  (3.15)  and  the  pressure 
derivative  in  x  to  equation  (3.16).  This  has  the  advantage  that  the  base 
pressure  distribution  is  fixed  and  disturbances  can  travel  out  of  the  domain. 
Similar  conditions  are  applied  in  this  work.  All  quantities  are  divided  into  a 
baseflow  and  disturbance  flow  at  the  free-stream 


<f>  =  he  +  <f>' .  (3.40) 

For  the  disturbances,  the  decay  condition  is  applied  and  the  base  part  is 
kept  constant.  A  detailed. derivation  of  the  decay  condition  can  be  found 


in  Thumm  (1991).  Here,  only  the  final  equation  is  presented.  For  the  most 
general  case,  a  compressible  flow,  the  decay  condition  has  the  form 


In  the  limit 


=  -\/<*2(l  -  M&(1  -  c)2)  +  7!  <fr' . 


to  ss 


dy  ’ 


(3.41) 


(3.42) 


equation  (3.41)  reduces  to  the  incompressible  version  (c.f.  Kloker,  1993) 

^-  =  -V/a2+72  <j>'.  (3.43) 


Figure  3.4  Validation  of  the  decay-condition.  Amplitude  distribution  of  a 
TS-wave  versus  Resi.  BIA  and  BA:  calculations  by  von  Terzi  (2004)  with 
decay  condition  on  total  flow  variables;  IVALDC:  decay  condition  only  on 
disturbance  variables. 

Figure  3.4  shows  a  test  calculation  (see  appendix  C  for  the  computational 
parameters)  where  the  previously  described  decay  condition  is  employed.  As 
a  test  case,  the  primary  instability  behavior  of  an  incompressible  flow  oyer  a 
flat  plate  is  simulated.  The  amplitude  of  the  2D  disturbances,  that  are  in¬ 
troduced  into  the  computational  domain,  is  small,  so  that  a  T S  wave  in  the 
linear  stage  is  excited.  The  graph  illustrates  the  logarithm  of  the  downstream 
amplitude  distribution  of  the  u- velocity  normalized  with  its  minimum  Sim¬ 
ulations  by  von  Terzi  (2004)  are  utilized  as  comparison.  Case  BIA  agrees 
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very  well  with  Fasel  &  Konzelmann  (1990)  and  is  taken  as  the  result  which 
reproduces  the  correct  amplitude  distribution.  It  was  obtained  from  an  in¬ 
compressible  code.  For  higher  Reynolds  numbers  (based  on  the  displacement 
thickness),  the  coarse  and  fine  version  of  BA  vary  from  BIA.  Both  versions 
of  BA  are  conducted  with  the  compressible  code  which  is  also  used  for  the 
present  research,  but  with  a  different  implementation  of  the  decay  condi¬ 
tion  which  did  not  fix  the  pressure  gradient  at  the  free-stream  boundary. 
IVALDC ,  the  case  with  the  new  decay  condition  does  not  deteriorate  the 
solution  of  the  case  without  APG. 


Compressible  Flow 

Figure  3.2  illustrates  the  steps  which  have  to  be  taken  to  start  a  forced 
flow  simulation  with  APG.  First,  a  baseflow  without  any  pressure  gradient 
is  generated.  Then,  at  the  free-stream,  a  specified  pressure  distribution  is 
introduced  and  the  solution  converged  out.  Finally,  this  baseflow  serves  as 
initial  condition  for  the  unsteady  simulation.  These  three  steps  need  dif¬ 
ferent  free-stream  boundary  conditions  which  are  explained  in  the  following 
paragraphs. 

For  the  baseflow  calculation  without  APG,  a  characteristic  boundary  condi¬ 
tion  derived  by  Harris  (1997)  is  used 

dd> 

Ts  =  0-  P-44) 

All  variables  do  not  change  along  the  direction  of  the  characteristics  that  is 
denoted  by  c  and  given  by  the  Mach  angle. 


Figure  3.5  Adjustment  between  similarity  solution  and  DNS  solution.  The 
pressure  distribution  of  the  first  30  points  is  plotted. 

At  the  inflow,  the  adjustment  between  the  similarity  solution  and  the  DNS 


solution  causes  an  effect  similar  to  an  expansion  fan  which  leaves  the  domain 
along  characteristics  (figure  3.5).  If  the  pressure  is  fixed  to  prescribe  a  desired 
pressure  distribution  at  the  free-stream,  which  is  explained  later,  the  expan¬ 
sion  fan  is  reflected  and  makes  it  impossible  to  create  a  clean  baseflow  with 
APG.  Therefore,  the  grid  points  where  the  pressure  is  strongly  influenced 
by  the  expansion  fan  are  cut  off  the  computational  domain  and  a  smaller 
domain  is  used  for  the  generation  of  the  baseflow  with  APG. 


Figure  3.6  Pressure  distribution  at  the  free-stream  (pIC)  and  its  derivative 
(dpic/dx).  xz  is  the  location  where  the  linear  distribution  with  positive 
slope  starts  and  £4  —  transdx  is  its  end.  £4  +  transdx  can  serve  as  a  start 
location  for  a  linear  distribution  with  negative  slope  which  ends  at  £5.  Here 
no  distribution  with  negative  slope  is  used,  transdx  is  a  parameter  which 
can  be  specified  in  the  DNS  code  and  which  defines  the  interval  where  a 
quadratic  Bezier-eurve  is  employed. 

In  order  to  obtain  a  baseflow  with  APG,  the  pressure  has  to  be  fixed  at 
the  free-stream.  Figure  3.6  shows  one  example.  In  all  calculations,  a  con¬ 
stant  gradient  was  employed.  It  is  not  possible  to  prescribe  this  pressure 
distribution  throughout  the  domain  from  the  inflow  to  the  outflow,  because 
at  the  corner  between  inflow  and  free-stream  a  shock  would  be  generated 
which  would  travel  along  a  characteristic  to  the  wall.  Therefore  a  blending 
function  has  to  be  used  which  smoothly  connects  the  original  pressure  distri¬ 
bution  (from  the  calculation  with  zero  pressure  gradient)  with  the  new  APG 
distribution.  For  this  purpose  a  sin2-function  is  employed  that  starts  at  £0 
and  ends  at  £3.  The  solid  line  represents  the  pressure  along  the  £-direction 


at  the  upper  boundary  and  the  dashed  line  its  first  derivative.  From  xz  on, 
the  pressure  is  increased  until  xA  -  transdx  and  its  distribution  has  the  slope 
A_pi/(a:4 — x3).  transdx  is  a  parameter  of  the  DNS  code  which  specifies  half  of 
the  interval  which  is  used  to  connect  the  pressure  distribution  upstream  of  x4 
with  the  distribution  downstream  of  x4.  As  a  blending  function,  a  quadratic 
Bezier-curve  is  prescribed.  A  Bezier-curve  is  defined  by  a  set  of  control 
points  which  determines  its  shape.  It  is  a  polynomial  and  was  invented  by 
the  French  engineer  Pierre  Bezier  in  the  60ies.  The  domain  downstream  of 
x4  -  transdx  represents  a  buffer  domain  that  can  be  used  to  accelerate  the 
flow  by  decreasing  the  pressure  linearly  with  the  slope  -Ap2/(r5  -x4)  from 
location  x4  +  transdx  to  x5.  The  domain  for  the  forced  flow  calculation 
is  indicated  as  “domain  of  interest.”  After  the  baseflow  with  the  APG  is 
converged  out,  the  unusable  parts  of  the  domain  are  cut  off  (for  example  up¬ 
stream  of  £3  and  downstream  of  x4 -transdx ).  The  derivative  of  the  pressure 
shows  no  jumps  at  the  locations  where  the  blending  functions  connect  differ¬ 
ent  pressure  distributions  and  this  is  necessary  for  the  ^-momentum  equation 
of  the  governing  equations  (3.3)  since  there  the  derivative  is  calculated.  The 
other  flow  variables  are  treated  by  applying  equation  (3.44). 

The  unsteady  simulations  for  the  compressible  flow  are  carried  out  with  the 
same  free-stream  boundary  condition  as  for  the  incompressible  flow.  All 
variables  are  treated  by  employing  equation  (3.41). 

3.5  Buffer  Domains  at  Inflow  and  Outflow 

The  buffer  domain  technique,  especially  ramping  down  disturbances  before 
they  reach  the  outflow  or  even  the  inflow,  is  very  effective  to  avoid  reflections 
of  disturbance  waves  at  the  boundaries.  Kloker  (1993),  Eissler  (1995)  and 
Meitz  (1996)  used  the  ramping  method  in  their  investigations.  Here,  this 
technique  is  applied  mainly  for  the  outflow.  For  the  incompressible  valida¬ 
tion,  a  ramping  function  is  also  used  for  the  inflow.  In  the  buffer  domain 
all  quantities,  p,  pUi  and  Et  are  multiplied  by  a  function  which  ramps  the 
disturbance  values  to  zero,  such  that  the  total  quantities  reach  their  base 
(initial)  values: 

fr(x)  =  ftc(x)  +  c(0(m  -  fic(x )) .  (3.45) 

In  equation  (3.45),  f(x)  symbolizes  the  total  flow  quantities  and  fic  their 
base  values.  The  difference  of  f(x)  -  fic(x),  which  is  the  disturbance  value, 
is  multiplied  by  a  weighting  function  c(£).  fr(x)  is  the  final  value  in  the 
buffer  region.  Meitz  (1996)  applied  a  weighting  function  of  the  form 


where  £  is  defined  as 


for  the  inflow,  and 
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for  the  outflow.  xstart  and  xend  are  the  start  and  end  location  of  the  buffer 
domain. 

During  the  simulation,  the  weighting  function  is  applied  at  every  Runge- 
Kutta  timestep  so  that 


fne^teP(x)  =  ^  foldstep  {x)  (3.49) 


If  this  operation  is  performed  n  times,  then  a  wave  that  travels  through 
the  buffer  domain  with  the  flow  “experiences”  the  multiplication  with  c(£) 
n  times.  This  leads  to  an  “effective  ramping  function”  of  c(£)n.  The  time 
a  wave  requires  to  travel  through  the  buffer  domain  is  determined  by  its 
phase  speed  and  the  length  of  the  buffer  domain.  From  this  time,  n  can  be 
'  calculated.  If  n  becomes  large,  repeated  multiplications  of  c(£)  with  itself 
cause  an  increasing  slope  of  the  effective  ramping  function.  This  is  definitely 
not  desirable  because  a  too  steep  effective  ramping  function  can  again  lead 
to  reflections  of  disturbances.  For  all  calculations,  n  was  small  enough,  such 
that  the  ramping  function  defined  in  equation  (3.46)  could  be  applied  on 
buffer  domains  with  the  length  of  one  or  two  streamwise  wavelengths  Xx. 


3.6  Disturbance  Generation 

Harmonic  disturbances  are  introduced  to  the  boundary  layer  by  periodic 
blowing  and  suction  through  a  slot  in  the  plate.  The  slot  is  located  about 
one  Tollmien— Schlichting  wavelength  downstream  of  the  inflow  (see  figure 
3.1).  All  other  boundary  conditions  are  not  affected.  Note  that  the  pressure 
boundary  condition  of  equation  (3.33)  is  only  approximately  valid  for  a  non¬ 
zero  wall-normal  velocity,  however,  it  is  a  reasonable  approximation  as  long 
as  A(3  /Re  =  AF  C  1,  where  (3  is  the  forcing  frequency  and  A  the  forcing 
amplitude.  The  disturbance  velocity  v  is  given  as 

v(x,  t )  =  Avp(x)  sin  (fit) .  (3.50) 

For  this  study,  the  amplitude  is  ramped  up  from  zero  over  the  first  forcing 
period,  i.e.,  for  0  <  t  <  T  —  2  7r//3.  In  addition,  a  zero  net  mass  flux  enters 
the  computational  domain  due  to  the  prescribed  spatial  disturbance  profile 
vp(x)  (see  figure  3.7,  with  the  circles  indicating  the  values  for  a  16-points 


resolution  on  a  slot  of  width  Ax/2).  In  order  to  balance  considerations  with 
respect  to  receptivity  and  resolution  requirements,  the  slot  should  be  between 
Ax/4  and  Ax  long  with  a  minimum  of  eight  points  in  the  x-direction.  The 
slot  itself  should  be  far  enough  away  from  the  inflow  (see  appendix  B.1.4). 


Figure -3.7  Disturbance  profile  vp(x)  (reproduced  from  Harris,  1997). 


3.7  Numerical  Method  for  Solving  the  Governing  Equations 

A  spatial  DNS  model  is  used  with  a  fourth-order  Runge-Kutta  method  for 
time-advancement  and  “fourth-order”  split  finite  differences  in  the  x-  and 
y-directions.  In  the  ^-direction,  a  periodic  solution  is  assumed  and,  conse¬ 
quently,  a  Fourier  transformation  is  applied.  Note  that,  for  three-dimensional 
calculations,  the  nonlinear  terms  of  the  governing  equations  are  computed  in 
physical  space  and  are  then  transformed  back  to  spectral  space  for  differenti¬ 
ation  and  integration.  Variables  are  symmetric  over  one-half  of  the  spanwise 
wavelength  (except  for  w,  which  is  antisymmetric  over  this  distance).  Thus 
only  half  a  wavelength  in  the  ^-direction  needs  to  be  computed.  In  order  to 
avoid  aliasing  errors,  the  number  of  physical  planes  is  chosen  to  be  greater  to 
3/2  spectral  planes.  The  numerical  method  is  explained  in  detail  in  Harris 
(1997)  and  von  Terzi  (2004). 


4.  VALIDATION  CASES 


The  original  code,  written  by  Harris  (1997)  for  investigations  of  compressible, 
plane  wakes,  has  already  been  extensively  validated  (c.f.  von  Terzi,  2004). 
Here,  several  additional  validation  cases  demonstrate  that  the  code  is  able  to 
perform  stability  investigations  of  supersonic  flow  over  a  flat-plate  with  and 
without  adverse  pressure-gradient  (APG).  Nonetheless,  physical  and  com¬ 
putational  parameters  must  be  chosen  carefully  to  accurately  capture  the 
physics  of  the  flow  and  to  ensure  the  stability  of  the  numerical  method.  First 
two  cases  from  the  literature  are  recomputed.  The  references  (Eissler,  1995; 
Thumm,  1991)  performed  supersonic  stability  investigations-one  at  a  Mach 
number  lower  and  the  other  at  a  higher  one  than  for  the  present  study.  The 
validation  of  the  code  with  pressure  gradient  is  divided  into  three  sections. 
First,  the  code  is  validated  for  an  incompressible  flow.  For  comparison,  case 
A  of  Kloker  (1993)  is  recomputed.  In  the  following  section  the  LST  solver  of 
Mack  is  described,  especially  the  steps  that  were  required  to  get  a  solution 
for  a  flow  with  pressure  gradient.  Finally,  two  DNS  calculations  (one  without 
and  one  with  APG)  are  compared  with  the  corresponding  solutions  of  the 
linear  stability  solver  of  Mack. 


Figure  4.1  Linear  stability  diagram  for  Ma  =  1.6;  the  solid  horizontal  line 
denotes  case  Al  of  Thumm  (1991). 


4.1  Supersonic  Flow  at  Mach  1.6 


Case  A1  of  Thumm  (1991)  allows  for  validating  the  growth  of  a  three- 
dimensional  TS-wave  with  results  from  DNS  and  linear  stability  theory.  The 
integration  domain  ranges  from  Rx  —  200  to  Rx  =  980  with  a  forcing  fre¬ 
quency  of  (3  =  5.0025.  This  case  is  marked  in  the  corresponding  stability 
diagram  (see  figure  4.1).  The  computational  parameters  are  summarized  in 
appendix  C,  table  C.l.  As  can  be  seen  in  figure  4.2  the  amplitude  and  phase 
distributions  agree  with  the  results  of  Thumm  (1991).  The  amplification  rate 


Figure  4.2  Distributions  of  amplitude  (left)  and  phase  (right)  for  u',  T'  and 
p'  (from  top  to  bottom)  at  Rx  -  680,  for  case  Al  of  Thumm  (1991). 

-&i  is  plotted  in  figure  4.3.  It  matches  the  calculation  of  Thumm  (1991)  and 
LST  fairly  well  demonstrating  that  the  DNS-code  is  capable  of  simulating 
the  growth  of  instability  waves  in  a  supersonic  boundary  layer. 
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Figure  4.3  Comparison  of  amplification  rates  for  case  Al  of  Thumm  (1991). 
4.2  Supersonic  Flow  at  Mach  4.8 

Case  C  of  Eissler  (1995)  simulates  a  flow  under  atmospheric  conditions  with 
constant  wall-temperature.  The  computational  parameters  are  summarized 
in  appendix  C,  table  C.l.  To  keep  the  wall- temperature  constant,  the  wall 
has  to  be  “cooled,”  thus  stabilizing  the  first  mode  and,  at  the  same  time, 
destabilizing  the  higher  modes  (Mack,  1969).  Indeed,  the  second  mode  for 
the  chosen  frequency  of  F  =  10  x  10-5  is  the  most  unstable  as  can  be  seen 
in  the  stability  diagram  of  figure  4.4. 


Figure  4.4  Linear  stability  diagram  for  Ma  =  4.8,  the  dashed  straight  line 
denotes  case  C  of  Eissler  (1995). 


33 


Figure  4.5  Distributions  of  amplitude  (left)  and  phase  (right)  for  u'  and  T'\ 
data  from  Eissler  (1995, top)  and  current  investigation  (bottom). 

Eissler  (1995)  chose  this  case  to  show  that  computations  with  a  thermally 
perfect  gas  are  not  worth  the  additional  computational  cost  in  comparison 
to  computing  with  a  calorically  perfect  gas.  The  velocity  field  is  practically 
identical  while  the  temperature  distribution  of  the  computation  with  a  ther¬ 
mally  perfect  gas  ranges  slightly  above  the  one  with  the  calorically  perfect 
gas  (see  figure  4.5).  Since  thermodynamic  aspects  are  not  under  investiga¬ 
tion,  the  higher  computational  cost  calculating  with  a  thermally  perfect  gas 
is  not  justified.  For  this  reason,  all  further  computations  are  carried  out  with 
a  calorically  perfect  gas. 

Note  that,  with  respect  to  stability  considerations,  this  validation  case  is 
calculated  such  that  all  disturbances  in  the  streamwise  and  the  spanwise 
direction  have  the  same  frequency  (fundamental  resonance).  The  amplitude 
and  phase  distribution  of  the  velocity  and  temperature  disturbances  for  the 
2D  and  first  3D  mode  shown  in  figure  4.5  agree  well.  The  slightly  higher 
temperature  amplitude  in  the  baseflow  disturbance,  T(1,0),  is  likely  due  to 
convergence  problems  of  Eissler  (1995). 

The  fact  that  the  amplitudes  employed  by  Eissler  (1995)  and  the  location  of 
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Figure  4.6  Comparison  of  the  downstream  development  of  the  maximum 
amplitudes  of  pu'  for  case  C  of  Eissler  (1995). 

the  disturbance  slot  are  unknown,  complicates  the  comparison  of  the  down¬ 
stream  development  of  the  maximum  amplitudes.  Nevertheless,  the  max¬ 
imum  pu'  amplitudes  agree  perfectly  after  the  disturbance  slot  for  the  2D 
disturbance.  The  3D  disturbances  show  similar  behavior,  i.e.,  the  wiggles 
look  alike  and  the  fundamental  resonance  occurs  at  the  same  downstream 
location  (see  figure  4.6). 

4.3  Incompressible  Flow  With  Pressure-Gradient 
4.3.1  Basefiow  Simulation 

As  mentioned  earlier,  a  Falkner-Skan  similarity  solution  serves  as  initial  con¬ 
dition  for  the  generation  of  the  basefiow.  This  initial  condition  has  a  Hartree 
parameter  of  /?#  =  —0.181  because,  during  the  convergence  of  the  basefiow, 
the  shape-factor  obtained  from  the  DNS  drifts  to  the  theoretical  value  asso¬ 
ciated  with  =  —0.18.  This  is  due  to  the  domain  height  used  for  the  DNS 
which  is  twice  as  high  as  the  domain  height  of  the  simulations  of  Kloker 
(1993).  The  pressure  distribution  of  the  initial  condition  is  kept  constant 
at  the  free-stream.  The  flow  properties  for  case  ICVALFR  are  specified  in 
appendix  C,  table  C.2.  The  computational  domain  starts  at  x  =  1.585  and 
ends  at  x  =  5.844.  The  converged  basefiow  and  the  Falkner-Skan  solution 
are  plotted  in  Figure  4.7  and  4.8. 

Figure  4.7  shows  the  u- velocity  over  77.  The  u- velocity  is  normalized  by  its 
free-stream  value  and  rj  is  calculated  employing  equation  (3.13).  Since  this 


Figure  4.7  Comparison  between  Falkner-Skan  solution  and  the  baseflow  ob¬ 
tained  from  DNS.  The  u-velocity  profile  is  plotted  over  the  similarity  variable 
77- 


Figure  4.8  Comparison  between  Falkner-Skan  solution  and  the  baseflow  ob¬ 
tained  from  DNS.  Shown  is  the  u- velocity  at  the  free-stream  Over  x. 


profile  is  a  similarity  solution  it  does  not  change  in  the  downstream  direction. 
In  Figure  4.8  the  distribution  of  the  edge  velocity  Ue  over  x  is  plotted. 
The  Falkner-Skan  values  are  calculated  using  equation  (3.15).  The  base- 
flow  matches  the  theoretical  values  of  the  Falkner-Skan  solution  except  for 
a  slight  shift  which  is  due  to  the  initial  condition.  However,  it  is  very  im¬ 
portant  to  compare  more  than  only  these  variables  to  find  out  whether  the 
same  baseflow  is  predicted  as  by  Kloker  (1993).  A  very  significant  quantity  is 
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the  shape-factor  Hu,  which  is  defined  as  the  ratio  between  the  displacement 
thickness  Si  and  the  momentum  thickness  62,  thus, 


The  displacement  thickness  of  an  incompressible  flow  is  given  by  the  equation 

6l  =Io  ^~W)dy'  (4-2) 

It  accounts  for  the  growth  of  the  boundary  layer  and  defines  the  displacement 
of  the  potential  flow  at  the  free-stream  by  the  value  8\.  The  equation  for  the 
momentum  thickness  has  the  form 

(4'3) 

Since  the  flow  near  the  wall  experiences  friction,  it  loses  momentum  which 
is  captured  by  the  momentum  thickness.  The  ratio  of  these  quantities  is  a 
constant  for  a  laminar  flow  which  possesses  a  similarity  characteristic.  The 
reason  for  this  fact  is  the  similarity  characteristic  of  the  flow.  The  shape- 


Figure  4.9  Shape-factor  of  the  Navier-Stokes  baseflow  and  the  theoretical 
value  of  the  similarity  solution  with  =  —0.18. 

factor  obtained  from  the  DNS  (figure  4.9)  matches  the  theoretical  value  for  a 
flow  with  Ph  =  —0.18  between  x  =  2.2  and  x  =  4.2.  The  disturbance  slot  for 
the  unsteady  simulation  with  controlled  forcing  is  located  between  x  =  2.082 
and  x  =  2.293.  For  this  reason,  the  results  of  the  unsteady  calculation  with 


37 


forcing  are  only  compared  to  the  results  of  Kloker  (1993)  downstream  of 
x  =  2.2.  Downstream  of  x  =  4.2,  the  2D  wave  and  the  3D  wave  start  to 
interact  with  each  other  nonlinearly.  The  resolution  of  the  simulation  is  too 
coarse  to  reproduce  the  correct  physics  for  this  stage.  For  this  reason,  it  is 
not  attempted  to  match  Kloker’s  curves  obtained  from  the  high  resolution 
simulation  (case  G)  downstream  of  x  =  4.2. 

Before  the  results  of  the  forced  flow  simulation  are  discussed,  attention  is 
being  devoted  to  the  Hartree  parameter  (3h .  The  initial  condition  has  a 
constant  Hartree  parameter.  It  is  of  interest  to  see  if  this  value  is  conserved 
in  the  DNS  calculation,  fig  can  be  computed  in  a  post-processing  step  with 
its  compressible  definition.  A  program  that  reads  in  the  data  of  the  converged 
baseflow  was  written  and  tested  with  the  previously  described  baseflow.  The 
definition  of  /?#■  is 


2£  dUejx) 
PH  Ue(x)  dl; 


(4.4) 


with 


£(x)  f  Ue(x)pe(x)p,e(x)dx . 
Jo 


(4.5) 


For  incompressible  flows  with  constant  viscosity,  p,e  and  pe  are  independent  of 
x.  Equation  (4.5)  is  part  of  the  Levy-Lee  transformation  (c.f.  Schlichting  & 
Gersten,  2000)  which  can  be  used  with  equation  (4.4)  to  derive  a  self-similar 
compressible  boundary  layer  with  streamwise  pressure  gradient. 

Since  the  inflow  of  the  computational  domain  (see  figure  3.1)  is  shifted  by  x0 
from  the  leading  edge  of  the  flat  plate,  £  has  to  be  calculated  with 


£(®)=/  ue(x)pe(x)pe(x)dx- f£0,  (4.6) 

Jxo 


where 

rx  o 

£o  =  j  ue (x)pe  (x)pe (x )dx . 

(4.7) 

In  most  cases, 

the  distribution  of  Ue(x),  pe{x)  and  pe(x) 

are  not  known 

upstream  of  the  computational  domain.  That  is  the  reason  why  equation 
(4.4)  is  easier  to  apply.  If  the  Hartree  parameter  is  given  at  the  inflow,  £0 

has  the  value 

_ 

rydUin  > 

(4.8) 

with 

dUf  _  dUjr  dx 

(4.9) 

~  dx  d£ 

and 

dx  1 

(4.10) 

•  di  ~  uinPfpf  • 

For  (3$  equal  to  zero,  dU™/dx  is  also  zero.  For  this  case,  equation  (4.7)  leads 
to 

4o  =  x0 .  (4.11) 

Figure  4.10  shows  the  result  of  the  previous  mentioned  post-processing  tool. 
The  Hartree  parameter  of  the  DNS  agrees  very  well  with  the  theoretical  value 
of  the  initial  condition. 
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Figure  4.10  Hartree  parameter  distribution  over  x  for  the  baseflow  obtained 
from  DNS. 


4.3.2  Forced  Flow  Simulation 

As  comparison,  a  calculation  is  chosen  in  which  Kloker  (1993)  simulated 
aligned  A- vortices  generated  due  to  a  strong  pressure  gradient.  This  physi¬ 
cal  process  which  creates  aligned  A- vortices  is  called  fundamental  resonance. 
Three  ToUmien-Schlichting  waves,  one  2D  and  two  3D,  interact  which  each 
other.  When  the  2D  ToUmien-Schlichting  wave  reaches  large  amplitudes,  it 
experiences  a  3D  deformation  which  is  caused  by  the  3D  waves.  According  to 
linear  theory,  2D  waves  are  more  amplified  than  3D  waves  for  subsonic  flows. 
However,  when  the  2D  wave  has  a  large  amplitude,  energy  is  being  trans¬ 
ferred  to  the  3D  waves  and  their  amplification  rate  changes  significantly.  The 
amplitude  of  the  2D  wave  saturates  and  the  3D  waves  can  reach  amplitudes 
of  the  order  of  the  amplitude  of  the  2D  wave. 

The  forced  flow  calculation  of  ICVALFR  is  conducted  with  almost  the  same 
computational  setup  as  Kloker  (1993)  used  for  his  simulation.  Disturbances 
with  the  angular  frequency  /?  =  10.8  are  introduced  into  the  domain  between 
x  =  2.082  and  x  =  2.293.  One  streamwise  wavelength  \x  contains  25  points 
and  the  domain  height  is  about  four  boundary  layer  thicknesses  high  (Kloker 
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used  only  two  boundary  layer  thicknesses).  The  computational  parameters 
are  summarized  in  appendix  C,  table  C.2.  First,  2D  calculations,  where  only 


Figure  4.11  Amplitude  distribution  of  the  u- velocity  versus  x  for  the  3D 
calculation  with  APG  according  to  Kloker  (1993):  2D  (left)  and  3D  (right). 

the  2D  Tollmien-Schlichting  waves  were  exited,  were  performed  to  determine 
the  correct  domain  height.  The  results  of  these  simulations  can  be  compared 
to  Kloker’s  3D  simulation  since  the  development  of  the  2D  wave  for  small 
amplitudes  behaves  according  to  linear  stability  theory. 

Figure  4.11  (left)  illustrates  the  influence  of  the  domain  height  on  the  am¬ 
plitude  distribution  of  the  2D  wave.  The  calculation  with  four  boundary 
layer  thicknesses  as  domain  height  can  reproduce  Kloker’s  data.  The  results 
of  the  3D  calculation  match  Kloker’s  data  as  well.  The  amplitude  distribu¬ 
tion  of  the  3D  Tollmien-Schlichting  wave  shown  in  figure  4.11  (right)  attains 
the  same  amplification  as  the  3D  wave  in  Kloker’s  simulation  for  the  linear 
stage.  At  x  =  4.0  the  onset  of  the  fundamental  resonance  can  be  observed. 
Therefore,  the  slope  of  the  amplitude  distribution  of  the  3D  wave  deviates 
from  the  slope  predicted  by  linear  stability  theory.  At  x  =  4.25,  ICVALFR 
starts  to  deviate  significantly  from  Kloker’s  values.  This  is  due  to  the  buffer 
domain  which  starts  at  x  =  4.5. 

4.4  Linear  Stability  Theory  With  Pressure  Gradient 

The  linear  stability  solver  of  Mack  is  a  very  powerful  tool  with  many  options 
to  compute  the  linear  stability  behavior  of  different  flow  types.  It  contains 
an  additional  tool  that  can  generate  the  similarity  solution  for  a  compressible 
flow  over  a  flat  plate  without  a  pressure  gradient.  The  output  of  this  tool  is 
then  read  into  the  solver  and  the  complex  wavenumber  for  a  specified  local 


Reynolds  number,  frequency  and  wave  angle  is  calculated  from  this  initial 
guess.  Although  the  similarity  solution  of  a  compressible  boundary  layer 
has  a  different  definition  of  the  similarity  variable  tj  than  the  incompressible 
version,  Mack  still  used  the  incompressible  definition  to  transform  the  linear 
stability  equations  into  the  77  coordinate  system  (c.f.  Thumm,  1991)  with 


V  =  lRx-  (4.12) 

We  demonstrated  before,  that  for  Mach  3,  this  procedure  is  a  very  good 
approximation. 

Since  the  global  Reynolds  number  of  Re  =  100, 000  is  fixed  in  the  stability 
solver,  Mack  scaled  the  complex  wavenumber  by  the  factor  Rx/Re  to  be 
independent  of  the  Reynolds  number 


or* 


amack  —  a  — 
'Re'  ~a'Re 


(4.13) 


For  comparison  with  the  DNS  data,  the  complex  wavenumber  has  to  be 
rescaled  by  the  Reynolds  number  which  is  used  in  the  DNS  simulation,  i.e.. 


tDNS  __  n. mack  .D./VS' 

T  T  RDNS  ’  a< 


=  aTack 


ReDNS 
R?NS  - 


(4.14) 


Flow  ] 

aroperties: 

Initial  values: 

Ma 

Too 

Re 

Wall 

3 

103.6iT 

1578102 

adiabatic 

Rx 

F 

Oii 

Oir 

280 

6.0  x  10-5 
3.490  x  10~4 
2.712  x  IQ"2 

Table  4.1  Physical  properties  of  the  experiments  (Brown,  2002)  and  the 
initial  values  for  the  2D  linear  stability  analysis. 


How  can  a  stability  analysis  with  pressure  gradient  be  performed  without 
having  any  corresponding  similarity  solution?  The  solution  for  that  problem 
can  be  found  using  the  DNS  data.  The  DNS  calculations  provide  a  converged 
baseflow  with  APG.  This  data  must  be  transformed  into  the  format  which  is 
read  into  the  stability  solver.  A  program  was  written  in  order  to  retrieve  the 
necessary  information  from  the  DNS  data  at  specified  streamwise  locations 
and  to  create  the  input  data  for  the  stability  solver.  Since,  for  this  case,  77 
(according  to  equation  4.12)  is  not  a  similarity  variable,  the  linear  stability 
analysis  can  only  be  performed  locally.  For  every  streamwise  location  a  new 


input  profile  must  be  generated.  As  initial  guess,  the  complex  wavenumber 
of  a  previous  stability  analysis  of  a  similarity  profile  with  the  same  physical 
properties  (except  APG)  is  used.  To  validate  this  procedure,  a  2D  stability 
analysis  with  the  physical  properties  of  the  experiments  by  Graziosi  &  Brown 
(2002)  and  without  APG  was  performed  (see  table  4.1).  The  investigated 
profiles  were  the  similarity  solution  of  Mack’s  similarity  solver  and  the  input 
profiles  created  from  the  corresponding  DNS  data  by  the  procedure  discussed 
above.  In  table  4.1  the  flow  properties,  disturbance  frequency,  start  location 
and  the  initial  eigenvalues  are  listed.  The  stability  analysis  started  at  Rx  = 
280  and  ended  at  Rx  =  880.  The  initial  eigenvalues  are  the  local  solution  for 
position  Rx  =  280. 


Figure  4.12  Eigenfunctions  of  u- velocity  (top  left),  v- velocity  (top  right) 
and  pressure  (bottom  left)  disturbances  at  Rx  =  280  and  amplification  rate 
o.i  over  Rx  (bottom  right)  from  stability  calculations  using  the  similarity 
solution  and  the  DNS  profiles. 

Figure  4.12  shows  the  results  obtained  using  this  stability  analysis.  The 
eigenfunctions  of  u,  v  and  p  at  position  Rx  =  280  and  t  the  amplification 
rate  a,  over  Rx  are  plotted.  Like  for  every  numerical  simulation,  the  optimal 
domain  height  of  the  integration  domain  and  the  largest  spatial  stepsize 


that  still  gives  the  correct  physical  result  had  to  be  determined.  Therefore, 
for  the  DNS  baseflow  profile,  which  was  read  into  the  stability  solver,  a 
resolution  and  domain  height  study  was  carried  out.  It  was  observed  that  for 
a  large  number  of  points  in  the  boundary  layer  in  the  wall-normal  direction 
(more  than  200  points)  Mack’s  stability  solver  could  not  produce  a  correct 
amplification  rate  anymore.  This  fact  marks  a  limit  for  all  DNS  based 
linear  stability  investigations  for  the  present  research  and  it  could  be  the 
explanation  for  the  difference  of  approximately  7.5%  between  the  growth 
rate  using  the  similarity  solution  as  input  profile  and  the  growth  rate  using 
the  DNS  data  as  an  input  shown  in  figure  4.12.  The  difference  of  around  7.5% 
of  the  growth  rate  is  determined  at  the  location  of  the  maximum  growth  rate 
using  the  similarity  solution  as  input  and  is  based  on  this  maximum  value. 
However,  the  difference  could  also  be  caused  by  the  slight  discrepancy  in  the 
profiles  of  the  similarity  solution  and  the  DNS  data. 

4.5  Compressible  Flow  With  Adverse  Pressure-Gradient 


Figure  4.13  Linear  stability  diagrams  for  Ma  =  3,  =  103.6AT,  adiabatic 

wall:  =  0°  (left)  and  %  ~  65°  (right).  The  region  outside  of  the  contours 

represents  the  region  of  damped  frequencies. 

In  this  section  two  DNS  calculations,  one  with  and  one  without  APG,  are 
presented.  The  flow  properties  are  listed  in  table  4.1.  A  more  detailed  table 
which  summarizes  the  computational  parameters  for  the  DNS  can  be  found 
in  appendix  C,  table  C.2.  Both  calculations  were  performed  with  small  dis¬ 
turbance  amplitudes,  so  that  they  can  be  compared  to  LST.  3D  disturbances 
were  excited  in  the  disturbance  slot.  This  choice  was  motivated  by  the  sta¬ 
bility  properties  of  the  flow.  The  2D  linear  stability  diagram  is  shown  in 
figure  4.13  (left).  Two  instability  modes  are  present,  a  viscous  (Tollmien- 
Schlichting  type)  and  a  second  (inviscid  or  Mack)  mode.  The  latter  is  more 
amplified.  It  is  very  difficult  to  excite  only  one  particular  instability  mode 


without  exciting  the  other  one.  Preliminary  2D  DNS  revealed  that  both  in¬ 
stability  modes  were  always  present  in  the  flow  field  and  that  the  Mack  mode, 
which  was  not  of  interest,  was  only  weakly  damped.  In  contrast,  the  3D  lin¬ 
ear  stability  diagram  (figure  4.13,  right)  contains  only  the  first  (TS-type) 
mode.  This  simplifies  the  post-processing  of  the  computational  data  since 
it  is  not  necessary  to  decompose  the  flow  field  into  the  different  instability 
modes.  The  disturbance  waves  are  introduced  at  a  wave  angle  of  T  ~  65°. 
This  angle  is  chosen  because  it  represents  the  angle  of  the  most  amplified 
wave  (Mack,  1984).  For  one  frequency  the  spanwise  wavenumber  is  a  con¬ 
stant  and  the  streamwise  wavenumber  is  a  function  of  x.  For  this  reason,  the 
wave  angle  varies  and  cannot  have  the  value  of  \fr  =  65°  in  the  whole  com¬ 
putational  domain.  The  spanwise  wavenumber  which  is  used  for  the  linear 
stability  analysis  is  defined  by  the  following  equation  (Eissler,  1995) 

7  =  dFRe,  (4.15) 

where  d  is  a  constant  that  defines  in  which  range  the  wave  angle  is  located. 
In  Mack’s  stability  solver  the  streamwise  wavenumber  a  is  normalized  by  the 
ratio  Rx/Re  (see  section  4.4).  This  scaling  is  also  applied  to  the  spanwise 
wavenumber.  Therefore,  equation  (4.16)  can  be  modified  to 

lmack  =  dFRx .  (4.16) 

Table  4.2  lists  the  values  of  d  for  both  linear  stability  analyses,  d  is  chosen  in 
a  way  such,  that  T  is  equal  to  65°  at  the  Rx  locations  where  the  amplitude 
distributions  of  the  DNS  and  the  eigenfunction  of  the  LST  are  compared  to 
each  other. 


No  pressure  gradient: 

Pressure  gradient: 

d 

3.0 

m 

3.5 

Table  4.2  Constant  d  for  the  two  different  simulations. 


4.5.1  Simulation  Without  Pressure  Gradient 

Figure  4.13  (right)  shows  the  stability  diagram  for  the  simulation  without  any 
pressure  gradient.  The  domain  starts  at  Rx  =  250  and  ends  at  Rx  =  1893. 
Disturbances  of  frequency  F  =  3  x  10~5  are  introduced  into  the  domain 
between  Rx  =  594  and  Rx  =  655.  The  buffer  domain  starts  at  Rx  =  1722. 
The  computational  domain  is  very  long  since  at  the  location  Rx  =  1400, 
where  the  DNS  is  compared  to  the  LST  data,  influences  of  the  disturbance 


slot  and  non-parallel  effects  are  required  to  be  weak.  The  resolution  in  the 
x-  and  y-directions  is  chosen  such  that  one  streamwise  wavelength  \x  is 
resolved  by  eleven  points  and  at  the  end  of  the  domain,  100  points  are  within 
the  boundary  layer  in  the  wall-normal  direction.  The  domain  height  is  two 
boundary  layer  thicknesses  at  the  outflow.  This  domain  height  is  also  used 
by  Thumm  (1991)  for  his  linear  stability  investigations.  • 


Figure  4.14  Amplitude  distribution  of  the  u-velocity  (top  left),  v-velocity 
(top  right)  and  pressure  (bottom  left)  at  Rx  =  1400  and  amplification  rate 
a.i  over  Rx  for  T  ~  65°  (bottom  right). 


In  figure  4.14,  the  amplitude  distributions  obtained  by  the  DNS  are  compared 
to  the  eigenfunctions  of  the  LST.  In  the  plot  for  the  u- velocity  (top  left),  it 
is  visible  that  for  a  3D  wave  the  eigenfunction  has  a  second  maximum  close 
to  the  wall.  Next  to  the  location  of  the  second  phase  shift,  the  amplitude 
distribution  from  the  DNS  varies  slightly  from  the  LST.  This  discrepancy 
was  also  observed  by  Thumm  (1991).  The  u-velocity  from  the  DNS  shows 
a  large  difference  to  the  LST  results  for  4  <  77  <  6.  The  deviation  becomes 
smaller  when  the  DNS  baseflow  is  used  for  the  linear  stability  analysis  (see 
section  4.5.2).  In  figure  4.14  (bottom  right)  the  amplification  rate  a*  over 
Rx  is  plotted.  It  is  calculated  from  the  maximum  value  of  the  w- velocity 


perturbations  in  y.  The  value  based  on  the  LST  deviates  by  27%  from  the 
DNS  solution  and  the  deviation  increases  to  infinity  for  afNS  equal  to  zero. 
This  means  that  the  DNS  predicts  a  higher  amplification  rate  than  LST. 
Thumm  (1991)  concluded  that  the  non-parallel  effects,  which  are  caused  by 
the  growth  of  the  boundary  layer,  are  responsible  for  this  behavior.  But  he 
did  not  present  any  proof  for  this  statement. 


Figure  4.15  Amplification  rate  a*  of  the  resolution  study  over  Rx  for  #  ~ 
60°  (left)  and  Fourier-transformed  u- velocity  in  time  (right,  the  amplitude 

distribution  between  points  100  - 150  in  x  and  points  1 - 80  in  y  is 

plotted). 

To  validate  if  the  DNS  values  are  correct  a  resolution  study  was  carried  out. 
The  results  are  shown  in  figure  4.15  (left).  All  cases  (with  different  resolution 
in  the  x-  and  ^-directions,  larger  domain  height,  smaller  disturbance  ampli¬ 
tude  and  different  inflow  buffer  domains)  reproduced  identical  results,  thus, 
indicating  a  well-resolved  simulation.  The  saw-toothed  shape  of  the  curves 
in  figure  4.15  is  an  artifact  of  the  post- processing  and  can  be  eliminated  by 
interpolating  the  locations  of  the  u- velocity  maxima.  Figure  4.15  (right)  il¬ 
lustrates  that  the  original  DNS  data  after  the  Fourier-transformation  in  time 
has  indeed  a  very  smooth  distribution. 

4.5.2  Simulation  With  Adverse  Pressure  Gradient 

Figure  4.16  shows  the  linear  stability  diagram  for  the  compressible  validation 
case  with  APG.  The  dashed-dotted  line  represents  the  neutral  curve  of  figure 
4.13  (right).  Only  the  region  between  Rx  =  1000  and  Rx  =  1600  includes 
an  APG.  For  the  case  with  APG,  the  neutral  curve  is  shifted  to  higher  fre¬ 
quencies  and  the  highest  amplification  rate  cxi  is  approximately  three  times 
larger  than  for  the  case  with  0H  =  0.  fa  is  not  a  constant  and  it  is  plotted 
in  figure  4.17.  The  value  of  fa  is  calculated  with  the  post-processing  tool 


V 


L«v«i 

* 

10 

OJGOOOOO 

9 

-0.000462 

8 

-0.000925 

7 

-0.001387 

6 

^001850 

•5 

HMJ02S12 

4 

•0.002774 

■9 

■0.003237 

2 

-0003699 

1 

•01004162 

Figure  4.16  Stability  diagram  for  Ma  =  3,  =  103.6AT,  adiabatic  wall  and 

4/  ~  65°.  Dashed- dotted  line  indicates  the  neutral  curve  without  APG.  The 
region  between  Rx  —  1000  and  Rx  =  1600  is  the  stability  diagram  with  APG. 


Figure  4.17  j3H  for  compressible  validation  case  with  APG  (case  CVALPG). 

described  in  section  4.3.1.  The  pressure  distribution  of  the  DNS  calculation 
is  plotted  in  figure  4.18  (left).  In  §  3.4.4,  it  was  already  discussed  how  a 
converged  baseflow  with  adverse  pressure-gradient  is  obtained.  A  pressure 
distribution  is  introduced  on  a  baseflow  through  the  upper  boundary.  This 
baseflow  is  converged  out  and  then  the  regions  which  are  not  needed  are  cut 
off.  Figure  4.18  (right)  shows  the  original  pressure  distribution. 

The  computational  domain  of  the  forced  flow  calculation  starts  at  Rx  —  899 
and  ends  at  Rx  =  1822.  Disturbance  waves  with  the  angular  frequency  of 
F  =  3  x  1CT5  are  introduced  through  the  disturbance  slot  located  between 
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Figure  4.18  Pressure  distribution  of  the  DNS  calculation  (case  CVALPG): 
after  (left)  and  before  (right)  the  computational  domain  is  cut;  only  every 
fourth  point  is  shown  for  clarity. 

Rx  —  1048  and  Rx  =  1083.  The  buffer  domain  starts  at  Rx  =  1640.  The 
resolution  in  x  and  y  is  the  same  as  for  the  calculation  without  APG  (case 
CVALNP).  Therefore,  only  ten  points  resolve  one  streamwise  wavelength  Xx 
and  80  points  in  the  wall-normal  direction  are  within  the  boundary  layer 
at  the  end  of  the  computational  domain.  .Compared  to  the  case  without 
pressure  gradient  (case  CVALNP),  the  boundary  layer  for  case  CVALPG  is 
thinner.  Incompressible  flow  behaves  very  differently.  An  APG  always  thick¬ 
ens  an  incompressible  boundary  layer  in  order  to  conserve  its  mass  flux.  For 
compressible  flows,  an  APG  causes  the  density  to  increase  and,  consequently, 
produces  a  larger  mass  flux  in  the  boundary  layer  without  the  need  to  thicken 
it.  If  the  increase  in  density  is  strong  enough,  the  boundary  layer  can  even 
become  thinner. 

In  figure  4.19,  the  amplitude  distributions  from  the  DNS  are  compared  to 
the  eigenfunctions  obtained  using  LST.  The  u- velocity  (figure  4.19,  top  left) 
matches  the  theoretical  values  pretty  well.  The  discrepancy  between  the 
DNS  and  the  LST  results  for  the  v-velocity  at  r)  «  5  is  much  smaller  than 
for  case  CVALNP.  This  corroborates  the  previous  statement  that  it  makes 
a  difference  whether  the  DNS  baseflow  or  the  similarity  solution  is  used  for 
the  LST  analysis.  For  higher  values  of  rj,  the  amplitude  distribution  of  the 
DNS  deviates  from  the  theory.  After  a  domain  height  study  it  was  concluded 
that  this  deviation  is  caused  by  the  upper  boundary  condition  of  the  stability 
solver.  With  increasing  domain  height  this  difference  decreases.  In  figure  4.19 
(bottom  right)  the  amplification  rate  a,  versus  Rx  is  plotted.  Like  for  case 
CVALNP,  the  waves  predicted  by  the  LST  calculations  are  amplified  about 
25%  less  than  the  waves  obtained  by  employing  DNS. 

Figure  4.20  shows  the  results  of  the  domain  height  study  for  the  stability 
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Figure  4.19  Amplitude  distribution  of  the  u-velocity  (top  left),  u-velocity 
(top  right)  and  pressure  (bottom  left)  at  Rx  =  1510  and  amplification  rate 
Oii  over  Rx  for  ^  ~  65°  (bottom  right). 


Ft, 


Figure  4.20  Amplification  rate  a{  versus  R?.  Domain  height  study  for  the 
hnear  stability  solver  of  Mack  (T  ~  61°). 


solver  of  Mack.  The  graph  indicates  that  the  domain  height  has  a  significant 
impact  on  the  amplification  rate  a*.  With  increasing  domain  height,  a, 
converges  to  the  value  that  is  obtained  when  the  domain  height  of  the  stability 
solver  is  chosen  to  be  two  boundary  layer  thicknesses. 


5.  PRINCETON  EXPERIMENTS 


At  Princeton  University,  Professor  Brown  and  co-workers  performed  an  ex¬ 
perimental  study  of  stability,  receptivity  and  transition  of  a  flat  plate  bound¬ 
ary  layer  at  Mach  3.  For  a  detailed  description  see  Graziosi  (1999)  and 
Graziosi  &  Brown  (2002).  In  the  following,  only  a  short  summary  of  their 
work  is  presented. 
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Figure  5.1  Quiet  region  of  LVTG  windtunnel  (reproduced  from  Graziosi  & 
Brown,  2002). 

The  Low  Turbulence  Variable  Geometry  (LTVG)  wind  tunnel  with  a  turbu¬ 
lence  rate  of  0.11%  to  0.39%  provides  a  28.5”  long  and  8”  wide  flat-plate  in 
the  test  section.  However,  the  quiet  region,  i.e.,  the  region  before  a  shock 
wave  is  reflected  at  a  tunnel  wall  and  hits  the  boundary  layer  again,  is  about 
9.5”  long  and,  consequently,  reaches  only  a  maximum  local  Reynolds  number 
of  Rx  =  725  (see  figure  5.1). 

Mean  flow  investigations  were  performed  at  a  stagnation  pressure  of 

ps  =  5 psia 

whereas,  due  to  the  lower  turbulence  level,  the  transition  investigation  was 
carried  out  at  a  stagnation  pressure  of 

ps  =  4 psia. 

For  a  summary  of  the  operating  conditions  see  table  5.1. 

Hot  wire  measurements  were  performed  at  five  different  downstream  locations 
between  Rx  =  400  and  Rx  =  800  along  the  centerline  of  the  plate.  The  data 
acquisition  locations  were  chosen  from  the  location  of  the  critical  Reynolds 
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Table  5.1  Operating  condition  of  LTVG  and  measuered  free-stream  turbu¬ 
lence  level. 

number  on  downstream.  Intensive  calibration  and  mean  flow  surveys  assured 
a  relatively  small  error  in  the  measured  data. 

Natural  disturbances  coming  from  the  upper  edges  of  the  wind  tunnel  pen¬ 
etrated  the  boundary  layer  and  excited  disturbances.  For  recent  studies  a 
loudspeaker  was  located  in  the  free-stream  to  introduce  controlled  distur¬ 
bances.  For  both  cases,  free-stream  turbulence  was  present  influencing  the 
measurements.  Frequencies  under  investigation  were: 

•  F  =  1.4  x  10-5-begin  of  the  instability  region  (branch  I), 

•  F  =  5  x  10_5-region  of  maximum  amplification  of  disturbances, 

•  F  =  8.1  x  10_5-end  of  the  instability  region  (branch  II), 

i.e.,  3,000  Hz ,  10,500  Hz  and  17,000  Hz,  respectively.  All  cases  under  inves¬ 
tigation  are  illustrated  in  the  linear  stability  diagram  of  figure  5.2. 


Figure  5.2  LST-diagram  with  investigated  cases  (reproduced  from  Graziosi, 
1999). 


6.  RESULTS  AND  DISCUSSION 


6.1  Baseflow  Simulations 

All  computations  are  carried  out  on  a  grid  according  to  the  results  of  the 
grid,  convergence  study  in  appendix  B.l,  i.e.,  the  domain  length  is  twelve  \x 
of  the  corresponding  TS-wave,  the  domain  height  is  nine  S^tfiow  with  an 
equidistant  grid  for  the  first  two  boundary  layer  thicknesses  and  a  stretched 
grid  above.  The  disturbance  slot  is  located  one  wavelength  downstream  of 
the  inflow  boundary  and  is  1/2  Xx  long  (exemptions  are  clearly  stated).  The 
resolution  is  32  points  per  Xx  and  ten  points  per  <5in/iou,. 


6.1.1  Comparison  With  Similarity  Solution 

The  difference  between  the' similarity  solution  and  the’ baseflow,  i.e.,  the 
steady  solution  to  the  Navier-Stokes,  equations  should  be  small.  This  con¬ 
jecture  is  corroborated  in  the  following.  First,  the  boundary  layer  thickness 
is  compared  in  table  6.1,  where  6*  denotes  the  dimensional  boundary  layer 
thickness.  As  reference  case  serves  the  boundary  layer  thickness  after  von 
Mises  (c.f.  Schlichting  &  Gersten,  2000)  for  a  laminar,  zero  pressure-gradient, 
incompressible  flat-plate  boundary  layer: 


S*  =  5 


(6.1) 


5*  [10-3m] 

von  Mises 

Navier-Stokes 

similarity  solution 

i?*  =467 

1.0711 

1.99097 

1.962 

i?*  =636 

1.4872 

2.68348 

2.6546 

A*  =769 

1.76371 

3.34714 

3.2029 

1! 

00 

4^ 

1.981655 

3.6645 

3.6357 

Table  6.1  Boundary  layer  thickness  comparison. 

Table  6.1  shows  that  the  estimate  according  to  von  Mises  is  only  a  rough 
guess  for  compressible  flows,  while  the  small  difference  between  the  similarity 
solution  and  the  Navier-Stokes  solution  is  most  likely  due  to  the  presence 
of  a  small  favorable  pressure-gradient  (see  below) .  The  maximum  deviation 
within  the  flow  field  is  about  3%  in  the  u- velocity.  In  figure  6.1,  the  similarity 
solution  and  the  baseflow,  normalized  with  the  value  at  the  edge  of  the 
boundary  layer,  are  compared.  The  increased  boundary  layer  thickness,  as 
already  shown  in  table  6.1,  can  be  observed  for  the  u- velocity.  While  the 


Figure  6.1  Comparison  of  Navier-Stokes  and  similarity  solution:  ^-velocity 
(top  left),  temperature  (top  right),  density  (bottom  left)  and  the  pressure  dis¬ 
tribution  of  the  Navier-Stokes  solution  in  the  streamwise  direction  (bottom 
right). 

density  p  shows  only  a  minor  deviation,  the  temperature  at  the  wall  is  larger 
for  the  Navier-Stokes  baseflow  than  for  the  similarity  solution.  However,  the 
difference  at  the  wall  is  only  6 K.  This  discrepancy  is  considered  to  have  an 
insignificant  impact  on  the  stability  characteristics  of  the  flow.  Nevertheless, 
the  difference  in  the  density  and  temperature  results  in  a  favorable  pressure 
gradient  in  the  streamwise  direction  because  the  equation  of  state  governing 
the  pressure  still  has  to  be  valid.  Figure  6.1  (bottom  right)  shows  a  typical 
pressure  distribution  in  streamwise  direction  while  the  pressure  in  the  y- 
direction  is  constant. 

By  scaling  data  at  three  different  downstream  locations  using  the  similarity 
variable  rj,  figure  6.2  shows  that  all  three  profiles  perfectly  match.  This 
demonstrates  that  the  flow  is  indeed  self-similar.  Matching  the  pressure 
gradient  in  the  region  of  interest  with  a  Falkner-Skan  similarity  velocity 
profile  led  to  a  Hartree-parameter  of  4.2  x  10-3  which  is  close  enough  to  zero 
to  assume  a  zero  pressure-gradient  flat-plate  boundary  layer  flow.  However, 


Figure  6.2  Navier-Stokes  solution  profiles  in  similarity  variable  rj  for  different 
downstream  locations. 

using  profiles  of  the  the  Navier-Stokes  solution  as  a  baseflow  for  a  LST 
analysis  results  in  a  smaller  unstable  region  due  to  the  favorable  pressure 
gradient  (see  figure  6.3).  One  could,  therefore,  conclude  that  the  Navier- 
Stokes  baseflow  should  be  more  stable  with  respect  to  disturbances  than  the 
similarity  solution.  However,  a  later  discussion  will  show  that  this  statement 
is  not  generally  true  (see  §  6.2.4). 


Figure  6.3  LST  diagram  of  similarity  and  Navier-Stokes  solution  for  T  =  60°. 
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R=677 


1^=836  R,=941 


Figure  6.4  Comparison  of  computation  and  experiment  at  different  down¬ 
stream  locations. 

6.1.2  Comparison  with  Experiments 

In  the  experiments,  baseflow  investigation  at  a  stagnation  pressure  of  5  psia 
are  performed.  Hot  wire  measurements  of  the  mass  flux  at  different  down¬ 
stream  locations  and  the  computed  streamwise  momentum  are  compared  in 
figure  6.4.  The  very  good  agreement  also  confirms  that  the  baseflow  simu¬ 
lation  captures  the  mean  flow  properties  at  Mach  3  correctly.  The  slightly 
different  profiles  farther  downstream  stem  from  the  reflected  shock  wave  (see 
figure  5.1).  In  the  experiments,  the  pressure  distribution  (figure  6.5)  shows 
only  small  pressure  deviations  in  the  streamwise  direction,  but,  because  of 
the  uncertainty  in  the  experimental  measurements,  an  overall  (favorable) 


pressure-gradient  could  still  be  present,  explaining  the  good  agreement  with 
the  computations. 
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Figure  6.5  Pressure  distribution  along  the  middle  axis  of  the  plate  (repro¬ 
duced  from  Graziosi,  1999). 


6.2  Linear  Stability  Behavior 

6.2.1  Disturbance  Frequency  F  =  1.4  x  10-5 

From  a  computational  point  of  view,  the  disturbance  frequency  of  F  =  1.4  x 
10-5  (/*  =  3, 000Hz)  is  problematic  because  a  very  large  TS-wavelength 
of  Xx  =  0.127m  is  linked  to  this  set  up.  Because  of  the  space  needed  for 


Figure  6.6  LST  stability  diagram  for  \fr  =  60°  and  'P  =  0°. 
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disturbance  generation  (see  §  3.6),  the  inflow  boundary  approaches  the  virtual 
origin  (x0  — *•  0)  in  figure  3.1  and  a  very  fine  computational  grid  has  to  be 
used  in  order  to  maintain  numerical  stability.  To  reduce  the  computational 
cost,  computations  on  an  immediately  stretched  grid  were  carried  out  and 
the  forcing  slot  was  located  only  15  grid  points  from  the  inflow  boundary 
with  an  extremely  small  width  of  only  five  grid  points.  •  Even  using  this 
disturbance  generation  configuration,  there  is  only  half  a  wavelength  space 
in  the  downstream  direction  before  the  disturbance  wave  enters  the  unstable 
region  (branch  I).  For  all  other  computational  parameters,  see  appendix  C, 
table  C.3. 

Comparison  with  Experiments 

To  check,  if  the  right  disturbances  are  introduced,  the  TS-wavelength  ob¬ 
tained  in  the  computations  is  compared  to  the  experimental  data  in  table  6.2. 
The  result  is  off  by  a  single  Ax  indicating  that  the  correct  TS-wave  is  cap¬ 
tured. 


K 

computation  [m] 
experiments  [m] 

0.127 

0.122 

Table  6.2  Wavelength  comparison  of  computation  with  experiments  for  F  = 
1.4  x  10~5. 


Figure  6.7  Amplitude  distribution  comparison  of  experiment  with  computa¬ 
tion  for  F  =  1.4  x  10~s  (/*  =  3,  OOOife). 

In  figure  6.7,  the  amplitude  distribution  in  the  streamwise  direction  is  plot¬ 
ted.  The  amplitude  of  the  introduced  disturbance  initially  decreases  rapidly. 


Therefore,  the  flow  is  not  fully  developed  causing  the  wiggles  in  the  ampli¬ 
tude  distribution.  So  it  is  questionable  if  the  physically  right  amplitude  at 
Rx  =  440  is  used  for  normalization.  In  both  cases,  the  amplitude  grows  at 
about  Rx  =  500  which  is  earlier  than  predicted  by  LST.  Farther  downstream, 
the  computation  shows  a  smaller  amplitude  than  the  experiments,  probably 
caused  by  the  normalization.  Nevertheless,  the  data  match  reasonably  well. 


Figure  6.8  Maximum  ^'-amplitude  (left)  and  amplification  rate  (right)  over 
Rx  for  F  =  1.4  x  10"5. 


Comparison  with  LST 

The  region  of  interest  is  too  close  to  the  disturbance  slot  (see  figure  6.6). 
Due  to  the  reasons  mentioned  above,  this  fact  cannot  be  changed  easily.  As 
a  consequence,  the  TS-wave  is  still  developing  which  causes  the  waviness 
of  the  amplitude  and  amplification  rate  plots  in  figure  6.8.  Instead  of  the 
amplification  rate,  the  maximum  it'-amplitude  is  therefore  used  to  give  some 
indication  of  the  flow  instability.  The  maximum  amplitude  in  figure  6.8  shows 
a  growing  amplitude  from  Rx  =  500  which  coincides  with  the  experiments 
but  is  earlier  than  LST  (see  figure  6.6).  Due  to  the  lack  of  data  farther 
upstream,  it  is  inconclusive  what  causes  this  effect.  From  Rx  =  700  on,  the 
logarithmic  plot  (figure  6.8)  shows  linear  behavior,  i.e.  exponential  growth 
of  the  disturbance  amplitude,  which  is  consistent  with  LST. 

In  figure  6.9,  the  eigenfunctions  of  LST  do  not  differ  too  much-their  shape 
stays  the  same.  The  discrepancies  stem  mainly  from  the  different  instabil¬ 
ity  modes  present  during  wave  development.  The  phase  distributions  show 
an  earlier  phase  shift  for  the  similarity  solution  than  for  the  Navier-Stokes 
solution  because  of  the  thicker  boundary  layer  of  the  latter.  Comparing  the 
eigenfunctions  with  the  amplitude  distribution  from  the  DNS  computation, 
the  second  maximum  at  -q  =  10  is  considerably  larger,  while  for  p'  and  T\  the 
first  maximum  is  smaller.  The  shape  of  the  eigenfunctions  of  the  computa- 


tion  is  not  similax  to  the  two  LST-eigenfunctions,  because  at  Rx  =  700,  the 
disturbances  have  only  travelled  one  wavelength  downstream  so  that  the  first 
maximum  has  not  yet  fully  developed.  The  absolute  value  of  the  phase  distri¬ 
bution  is  different  but  the  phase  shifts  are  consistent  with  the  eigenfunctions. 


Figure  6.9  Comparison  of  amplitude  (left)  and  phase  (right)  distribution  of 
the  computation  to  LST  based  on  the  similarity  and  the  Navier-Stokes  .  solu¬ 
tion;  ahown  are  u- velocity  (top),  density  (center)  and  temperature  (bottom); 
F  =  1.4  x  10-5,  Rx  =  700. 


6.2.2  Disturbance  Frequency  F  =  5.0  x  10~5 

This  case  was  chosen  due  to  its  proximity  to  the  region  of  maximum  ampli¬ 
fication  rate.  For  computational  parameters,  see  appendix  C,  table  C.4. 

Comparison  with  Experiments 

Comparison  of  the  TS-wavelengths  gives  confidence  that  the  same  distur¬ 
bances  as  in  the  experiments  are  introduced  (see  table  6.3). 


computation  [m]  0.0376 

experiments  [m]  0.0376 


Table  6.3  Wavelength  comparison  of  computation  with  experiments  for  F  = 
5  x  lO"5. 

In  figure  6.10,  the  DNS  results  are  compared  with  the  experimental  data. 
Shown  are  the  amplitude  distribution  for  different  downstream  locations  and 
the  amplification  rate  at  Rx  =  500  for  different  frequencies.  The  simulations 
and  the  experiments  match  fairly  well  at  lower  Rx.  The  larger  amplitude  for 
the  computation  farther  downstream  suggests  that  non-linearities  are  present 
in  the  experiments-probably  because  of  the  reflected  shockwave  at  Rx  ~  700 
(see  §  5). 


Figure  6.10  Comparison  of  experiment  with  computation  for  F  =  5  x  10-5: 
amplitude  distribution  over  Rx  (left)  and  amplification  rate  at  Rx  =  500  for 
different  frequencies. 


-0.0005 


Figure  6.11  Amplification  rate  comparison  of  computations  with  results  from 
LST  for  F  =  5  x  10"5  (/*  =  10, 500ifc). 

Comparison  with  LST 

The  computations  and  the  experimental  data  both  show  a  larger  amplifica¬ 
tion  rate  than  LST,  where  the  amplification  rate  of  the  computations  is  even 
larger  than  in  the  experiments  (see  figure  6.10,  right).  In  figure  6.11,  the 
maximum  amplification  rate  is  significantly  larger  compared  to  LST.  Farther 
downstream,  the  amplification  rate  curve  approaches  the  LST-curve  which 
can  also  be  seen  in  figure  6.11.  As  shown  in  figure  6.12  the  amplitude  distri¬ 
butions  of  the  computation  are  wider  compared  to  the  LST-eigenfunctions, 
both  with  similarity  and  Navier-Stokes  baseflows.  The  first  maximum  of  the 
u'-amplitude  distribution  of  the  computation  is  not  as  large,  but  all  eigen¬ 
functions  match  in  the  free-stream  ( r]  >  10).  The  phase  distributions  are 
similar  with  the  phase  shifts  consistent  to  the  eigenfunctions. 

6.2.3  Disturbance  Frequency  F  =  8.1  x  1CT5 

In  this  case,  the  branch  II  of  the  stability  diagram  for  a  wave  angle  T  =  60° 
is  crossed  at  Rx  =  700  with  F  =  8.1  x  10-5  (see  figure  6.6).  For  a  complete 
list  of  computational  parameters,  see  appendix  C. 


K 

computation  [m] 
experiments  [m] 

0.0241 

0.0249 

Table  6.4  Wavelength  comparison  of  computation  with  experiments  at  F  = 
8.1  x  10~5. 
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Figure  6.12  Comparison  of  amplitude  (left)  and  phase  (right)  distribution  of 
the  computation  to  LST  based  on  the  similarity  and  the  Navier-Stokes  solu¬ 
tion;  ahown  are  u- velocity  (top),  density  (center)  and  temperature  (bottom); 
F  =  5  x  10"5,  =  700. 


Comparison  with  Experiments 

As  above,  the  TS-wavelength  (table  6.4)  is  used  to  confirm  the  generation 
of  the  same  kind  of  disturbances  as  in  the  experiments.  The  good  agreement 
verifies  the  capturing  of  the  right  disturbance  frequency  and  assures  that  no 
sound-wave  or  numerical  noise  is  amplified  to  the  order  of  the  wave  under 
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Figure  6.13  Amplitude  distribution  comparison  of  experiment  with  compu¬ 
tation  for  F  =  8.1  x  10~5  (/*  =  17,  000Hz). 

investigation.  The  amplification  rate  in  figure  6.10  (right)  for  the  compu¬ 
tations  is  about  10%  larger  than  Graziosi  (1999)  found  in  his  experiments. 
Taking  the  uncertainty  of  supersonic  experiments  into  account,  the  compu¬ 
tations  match  the  experiments  reasonably  well.  In  figure  6.13,  the  amplitude 
distribution  in  downstream  direction  is  shown.  The  overall  agreement  is  re¬ 
markable.  Farther  downstream,  the  computation  shows  a  higher  amplitude 
than  the  experiments.  This  could  again  be  due  to  nonlinear  interaction  which 
is  avoided  in  the  computation. 


Figure  6.14  Amplification  rate  comparison  of  u'  and  pu'  with  LST  for  F  = 
8.1  x  10-5  (/*  =  17,  OOOife). 
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Figure  6.15  Comparison  of  amplitude  (left)  and  phase  (right)  distribution  of 
the  computation  to  LST  based  on  the  similarity  and  the  Navier-Stokes  solu¬ 
tion;  ahown  are  u- velocity  (top),  density  (center)  and  temperature  (bottom); 
F  =  8.1  x  1(T5,  Rx  =  700. 

Comparison  with  LST 

Although  only  small  disturbances  of  v'  are  introduced,  the  larger  amplifi¬ 
cation  rate  and  the  shift  of  the  instability  region  farther  downstream  are 
significant  (see  figure  6.14).  For  detailed  discussion  see  §  6.2.4.  Looking 
at  the  amplitude  distribution  in  figure  6.15,  the  LST  based  on  the  similar- 
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ity  solution  overestimates  the  temperature  and  density  amplitude,  while  the 
eigenfunctions  of  the  computations  are  closer  to  the  eigenfunctions  computed 
with  the  Navier-Stokes  solution.  Because  the  differences  between  the  base- 
flows  obtained  using  the  similarity  and  Navier-Stokes  solutions  are  small 
(see  §  6.1.1),  it  is  unclear  what  causes  the  large  deviation  in  the  p'-  and  T'- 
eigenfunction  of  the  LST  similarity  solution.  The  phase  distribution  of  the 
computation  shows  the  same  properties  than  the  LST-phase  distributions, 
for  both  the  similarity  and  the  Navier-Stokes  baseflow. 

6.2.4  Non-Parallel  Effects 

The  growth  of  the  boundary-layer  in  the  downstream  direction  causes  a 
change  in  the  wavenumbers  of  the  downstream  traveling  TS-waves.  This 
effect  is  not  included  in  the  investigations  employing  linear  stability  theory. 
For  LST,  a  “parallel  flow  assumption”  is  used  which  inherently  assumes  that 
the  baseflow  cannot  change  in  the  downstream  direction.  In  the  following, 
we  will  investigate  how  the  growth  of  the  boundary  layer  changes  the  region 
of  instability,  i.e.,  whether  it  moves  the  start  (branch  I  of  the  neutral  curve 
in  a  linear  stability  diagram)  and  the  end  (branch  II)  of  the  unstable  region. 

6.2.4. 1  Investigation  of  Branch  I 

In  §  6.2.1,  the  simulation  with  frequency  F  =  1.4  x  10-5  was  carried  out 
such  that  the  generated  TS-wave  crosses  branch  I.  However,  due  to  the  fact 
that  the  space  between  the  end  of  the  disturbance  slot  and  the  beginning 
of  the  unstable  region  is  only  \x/2  in  downstream  direction,  the  TS-wave 
was  not  fully  developed  when  entering  the  unstable  region.  Nevertheless, 
the  amplitude  comparison  of  figure  6.8  suggests  that  an  increased  unstable 
region  for  branch  I  may  exist  due  to  non-parallel  effects  of  the  compressible 
boundary  layer.  This  is  corroborated  by  the  experiments  where  the  amplitude 
growth  is  significant  enough  to  exceed  any  error  in  the  measurements. 

6.2.4.2  Investigation  of  Branch  II 

To  determine  the  increase  of  the  unstable  region  near  branch  II  due  to  non¬ 
parallel  effects,  four  frequencies  are  considered: 

•  F  =  5  x  10-5,  from  R  =  700  to  R  =  1450 

•  F  =  8.1  x  10~5,  from  R  =  400  to  R  =  1000 

•  F  =  1.5  x  10~4,  from  Rx  =  100  to  R  =  650 

•  F  =  2  x  10-4,  from  R  =  100  to  R  =  550 


For  a  detailed  description  of  the  two  lower  frequencies,  see  §  6.2.2  and  §  6.2.3, 
where  the  lowest  frequency  is  shifted  downstream  to  cross  the  neutral  curve 
at  branch  II.  The  two  higher  frequencies  are  chosen  to  investigate  the  non¬ 
parallel  effects  on  high-frequency  disturbances.  Detailed  listings  of  used  pa¬ 
rameters  can  be  found  in  appendix  C. 

In  the  incompressible  case,  the  best  criterion  for  the  unstable  region  of  the 
stability  diagram  is  the  inner  maximum  of  u'  (Fasel  &  Konzelmann,  1990). 
Here,  the  second  maximum  of  a  given  parameter  is  chosen,  because  the  sec¬ 
ond  maximum  is  much  larger  than  the  first  maximum  for  the  compressible 
eigenfunction.  In  this  study,  different  criteria  are  presented  to  determine 
non-parallel  effects,  but  no  conclusion  is  drawn  on  how  non-parallel  effects 
are  measured,  i.e.,  how  the  downshift  of  the  neutral  point  is  related  to  non¬ 
parallel  effects.  Generally,  the  following  quantities  are  investigated: 

•  streamwise  velocity  v! 

•  wall-normal  velocity  v' 

•  density  p' 

•  temperature  T' 

•  pressure  p' 

•  spanwise  vorticity  u'z 

•  disturbance  mass  flux  in  streamwise  direction  pu' 

Independent  of  the  frequency,  the  parameters  u\  p\  T  and  pu'  are  all  ampli¬ 
fied  in  the  same  manner  while  v\  p'  and  to'  show  different  behavior.  To  ensure 
that  the  downshift  is  not  governed  by  the  influence  of  a  pressure-gradient  in 
the  simulation,  the  neutral  curve  of  the  LST  based  on  the  Navier-Stokes 
solution  is  also  computed.  The  Navier-Stokes  solution  has  a  smaller  unsta¬ 
ble  region  compared  to  the  LST  governed  by  the  similarity  solution  which 
is  consistent  with  a  favorable  pressure-gradient.  This  results  in  an  even  far¬ 
ther  downshift  of  the  neutral  point  which  seems  to  be  accurate  since  the 
computation  captures  the  experimental  results  very  well. 

The  long- wavelength  oscillations  in  p'  indicate  that  sound  waves  overlay  the 
amplification  rate.  These  oscillations  also  influence  pu'  and  T'.  The  small  os¬ 
cillations  are  introduced  by  the  post-processing.  Therefore,  the  amplification 
rates  are  averaged  to  obtain  the  neutral  point. 


Frequency  F  =  5  x  10~5 


The  lowest  frequency  is  closest  to  the  neutral  curve  computed  with  LST 
leading  to  the  conclusion  that  the  non-parallel  effects  are  less  pronounced. 
Looking  at  it',  p',  V  and  pu'  the  neutral  point  in  figure  6.16  (left)  is  shifted 
from  Rx  =  1220  to  Rx  =  1270,  that  is  a  difference  of  A Rx  =  50  or  1.5  A* 
compared  to  the  LST  results  employing  the  similarity  solution  as  baseflow. 
Compared  to  the  LST  using  the  Navier-Stokes  baseflow,  the  shift  is  from 
Rx  =  1025  to  Rx  =  1270  or  6.8  \x.  The  ^'-amplification  rate  and  the 
LST  similarity  solution  are  the  same  so  that  non-parallel  effects  could  not 
be  determined,  but  the  downshift  compared  to  the  LST  with  Navier-Stokes 
baseflow  is  from  Rx  =  1025  to  Rx  =  1220  (5.3  Xx).  The  ^'-amplification 
rate  also  takes  much  longer  before  it  is  fully  developed.  Even  worse  is  the 
p'-amplification  rate  which  gives  the  smallest  difference  to  the  LST  with 
Navier-Stokes  baseflow,  i.e.,  a  downshift  to  Rx  —  1140  or  a  difference  of 
three  A*  (this  point  may  still  vary  because  of  the  late  development). 


Figure  6.16  Amplification  rate  for  different  reference  quantities;  F  =  5  x  10-5 
(left)  and  F  =  8.1  x  10-5  (right). 


Frequency  F  =  8.1  x  10~5 

For  the  frequency  of  /*  =  17,000 Hz  (see  figure  6.16,  right ),  the  instability 
region  is  shifted  about  A Rx  =  105  (from  R^  =  700  to  R^  =  805)  or  3 
Ax  downstream  if  the  ^'-amplification  rate  and  equivalent  parameters  are 
compared  to  the  LST  similarity  solution.  Compared  to  the  LST  Navier- 
Stokes  solution,  the  shift  is  from  Rx  =  600  to  Rx  =  805  (5.3  A*).  Looking  at 
the  u'- amplification  rate,  the  shift  is  less  than  for  the  latter  parameters,  i.e. 
the  shift  is  about  A R^  —  50  and  ARX  =  150  compared  to  the  LST  similarity 
solution  and  to  the  LST  Navier-Stokes  solution  respectively.  A  relationship 
between  the  development  of  the  v'  and  the  u'z- -amplification  rate  is  clearly 


visible-once  the  ^'-amplification  has  developed,  the  vorticity  amplification 
rate  is  getting  smaller  until  it  is  parallel  to  the  other  amplification  rate  curves. 
The  same  mechanism  can  be  seen  for  all  frequencies  under  investigation, 
although  it  is  not  this  obvious.  For  the  frequency  of  F  =  8.1  x  10~5,  the  ex¬ 
emplification  stays  also  different  to  the  ^'-amplification  rate  curve  while  for 
the  other  frequencies  it  approaches  the  same  values  due  to  this  mechanism. 
The  pressure  is  again  less  amplified  than  the  LST  similarity  solution  predicts 
and  has  therefore  the  smallest  downshift  to  Rx  =  655  (1.5  Xx)  compared  to 
the  LST  Navier-Stokes  solution. 

Frequency  F  =  1.5  x  10~4 

Unfortunately,  the  second  highest  frequency  shows  some  slight  oscillations 
around  the  neutral  point  so  that  no  precise  location  can  be  specified  (see 
figure  6.17,  left).  Figure  6.17  (left)  also  shows  that  these  oscillations  are  not 
produced  by  a  reflected  expansion  fan  as  experienced  for  the  grid  conver¬ 
gence  study  in  appendix  B.1.1,  since  the  oscillations  appear  with  the  same 
severity  at  the  same  location  independent  of  the  domain  height.  It  is  more 
likely  that  they  are  caused  by  a  numerical  error  when  the  amplitude  is  close 
to  zero  during  the  post-processing  procedure.  Computations  at  a  lower  fre¬ 
quency  do  not  show  this  behavior  because  of  the  higher  amplification  rate. 
Nonetheless,  the  amplitude  distribution  over  the  flow  field  (see  figure  6.17, 


Figure  6.17  Amplification  rates  for  different  domain  heights  (left)  and  am¬ 
plitude  distribution  of  u'  (right);  F  =  1.5x  10~4  (/*  =  31,250 Hz). 


right)  is  smooth  so  that  the  physical  properties  of  the  flow  are  captured  cor¬ 
rectly.  The  comparison  to  LST  is  done  with  a  domain  height  of  20  boundary 
layer  thicknesses.  Considering  n'  and  the  equivalent  parameters  as  criterion 
to  determine  non- parallel  effects,  the  region  of  the  neutral  point  in  figure  6.18 
(left)  lies  therefore  between  Rx  =  465  and  Rx  =  500  revealing  an  enormous 


shift  downstream  compared  to  Rx  =  350  from  the  LST  similarity  solution. 
The  maximum  amplitude  distribution  leads  to  the  conclusion  that  the  neu¬ 
tral  point  is  considered  to  be  at  Rx  =  480  (see  figure  6.19,  left).  The  overall 
shift  is  therefore  A Rx  =  130  or  3  Xx.  The  difference  between  the  LST  sim¬ 
ilarity  solution  and  LST  Navier-Stokes  solution  is  only  A Rx  =  40  so  that 
the  difference  between  the  ^'-amplification  rate  and  the  LST  Navier-Stokes 
solution  is  about  A Rx  =  170  (4  lx).  The  ^'-amplification  rate  is  closer  the 
^'-amplification  rate  curve  than  for  the  lower  frequencies  leading  to  the  con¬ 
clusion  that  the  non-parallel  effects  are  stronger  for  v'  than  for  the  smaller 
frequencies.  Considering  v'  as  criterion  to  determine  non-parallel  effects  re¬ 
sults  in  a  downshift  of  A Rx  =  90  (2  Xx)  or  AR^  =  130  (2.5  A*)  for  the  com¬ 
parison  with  the  LST  employing  the  similarity  solution  and  Navier-Stokes 
baseflow,  respectively.  The  pressure  amplification  rate  develops  too  late  to 
determine  the  neutral  point,  but  because  of  the  difference  to  the  other  ampli¬ 
fication  rates  farther  downstream,  it  is  considered  to  give  the  best  agreement 
with  the  LST  using  the  Navier-Stokes  baseflow. 


Figure  6.18  Amplification  rate  for  different  reference  quantities;  F  =  1.5  x 
10~4  (left)  and  F  =  2  x  10~4  (right). 


Frequency  F  =  2x  10~4 

Figure  6.18  (right)  shows  that  for  this  frequency  the  Navier-Stokes  solution  is 
already  stable  according  to  LST.  Looking  at  the  amplitude  distribution  of  u' 
in  figure  6.19  (right),  the  neutral  point  is  at  about  Rx  =  345.  The  resulting 
shift  of  the  neutral  point  downstream,  as  shown  in  figure  6.18  (right),  is 
about  Rx  ~  100  or  2.5  Xx  compared  to  the  LST  similarity  solution  where  the 
neutral  point  is  at  Rx  =  250.  Because  of  the  later  development  of  the  v'  and 
p'-amplification  rate  no  data  to  determine  the  neutral  point  can  be  extracted. 
Taking  the  ^'-amplification  rate  into  account  the  the  neutral  point  is  slightly 
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Figure  6.19  Maximum  amplitude  distribution  of  is';  F  —  1.5  x  10-4  (left) 
and  F  =  2  x  10-4  (right). 

farther  downstream,  namely  at  Rx  =  360.  • 

6. 2.4.3  Summary  of  Observed  Non- Parallel  Effects 

Non- parallel  effects  are  most  likely  present  for  the  flat  plate  boundary  layer 
at  Ma  =  3.  Because  the  experiments  experience  effects  that  agree  qualita¬ 
tively  with  the  computations,  it  is  unlikely  that  the  downshift  is  caused  by 
a  numerical  effect.  But  the  far  downstream  shift  of  the  neutral  point  for  is', 
p',  T'  and  pv!  is  not  only  governed  by  non-parallel  effects.  This  tremendous 
shift  may  be  explained  by  an  interaction  of  a  higher  mode  which  is  for  this 
Reynolds  and  Mach  number  not  yet  unstable  (otherwise  it  would  be  present 
in  the  LST-diagram).  Figure  6.6  shows  an  upcoming  second  mode  for  a  2D 
disturbance. 


Figure  6.20  Verification  of  linear  disturbance  assumption. 

Another  reason  could  be  the  resolution  in  the  spanwise  direction  so  that, 
although  linear  disturbances  are  introduced,  the  energy  wants  to  shift  to 


higher  modes.  To  check  that  the  introduced  perturbations  are  small  enough 
and,  consequently,  the  linearity  assumption  holds,  the  already  small  forcing 
amplitude  is  reduced  by  a  factor  of  100  (see  figure  6.20).  The  perfect  align¬ 
ment  of  the  curves  demonstrates  that  the  disturbances  behave  linearly  and 
the  comparison  with  LST  is  valid.  Therefore,  the  spanwise  resolution  is  not 
considered  to  influence  the  stability  behavior. 

The  good  agreement  between  the  ^'-amplification  rate  and  the  LST  similarity 
solution  is  remarkable  leading  to  the  conclusion  that  the  other  mechanism 
influencing  the  stability  behavior,  whatever  that  mechanism  is,  has  less  in¬ 
fluence  on  the  v'- velocity.  Nonetheless,  the  downshift  compared  to  the  LST 
Navier-Stokes  solution  is  remarkably  large. 


Figure  6.21  Pressure  eigenfunction  at  Rx  =  940. 

Without  proof,  the  p'-amplification  rate  compared  to  the  LST  with  Navier- 
Stokes  baseflow  is  probably  best  suitable  to  determine  non-parallel  effects. 
Mack  categorizes  the  modes  by  the  zero  crossings  of  the  pressure  eigenfunc¬ 
tion  (Mack,  1984).  Because  no  zero  crossing  is  visible  in  the  computations 
(see  figure  6.21),  the  pressure  distribution  is  associated  with  the  first  mode, 
i.e.,  the  TS-mode.  Therefore,  it  is  less  likely  that  higher  modes  are  influencing 
the  amplification  of  the  pressure  disturbances. 


Figure  6.22  Boundary  layer  growth  and  development  of  p'- wave  over  x. 

This  study  shows  that  independent  of  the  criterion,  non-parallel  effects  are 
much  stronger  than  in  the  incompressible  case  (c.f.  Fasel  &  Konzelmann, 


1990).  To  demonstrate  that  the  result  is  still  plausible,  figure  6.22  shows  the 
boundary  layer  growth  over  the  ^-direction  in  comparison  to  the  development 
of  the  wall  pressure  disturbance  wave.  The  fast  growing  boundary  layer  and, 
consequently,  the  increasing  u-velocity  indicate  that  the  assumption  of  a 
locally  parallel  flow  used  in  LST  is  not  necessarily  fulfilled  for  the  current 
investigation. 

6.3  Nonlineax  Resonances 

6.3.1  Classical  Fundamental  Resonance 


Figure  6.23  Fundamental  Breakdown.  Amplitude  distribution  of  u- 
perturbation  versus  downstream  direction.  Re  =  105,  F  =  5  x  10~5, 
=  60°. 

Known  from  incompressible  investigations  the  fundamental  and  subharmonic 
resonances  are  both  strong  mechanisms  to  transition.  Now  for  the  supersonic 
boundary  layer,  our  simulations  should  shed  some  light  on  the  question  if 
they  are  also  a  relevant  mechanism  for  Ma  =  3.  Several  computations  were 
conducted  to  reproduce  a  fundamental  and  subharmonic  resonance.  These 
simulations  indicated  that  the  largest  secondary  growth  rate  occurs  between 
a  two-dimensional  disturbance  and  a  disturbance  with  the  wave  angle  close 
to  the  one  associated  with  the  linearly  most  amplified  wave. 

For  a  fundamental  resonance  to  take  place  within  the  computational  domain, 
the  “classical”  fundamental  breakdown  needs  a  disturbance  amplitude  of  10% 
of  the  free-stream  velocity  (see  figure  6.23).  Due  to  these  large  amplitudes, 
the  fundamental  breakdown  does  not  seem  to  play  an  important  role  for  Ma 
=  3  transition.  It  is  to  the  authors  opinion  that  the  high  amplitudes  levels  are 
necessary  because  of  the  small  linear  amplification  of  two-dimensional  first 
and  second  mode  waves  (c.f.  figure  6.6).  An  APG  only  weakly  influences 


Figure  6.24  Isosurface  of  Q- criterion  Q  =  2  (left)  and  spanwise  vorticity  ioz  = 
15  (right)  for  a  classical  fundamental  resonance;  Re  =  10s,  Rx  =  740  -  895, 
F  =  5  x  10"5,  =  60°. 


the  growth  rates  of  a  two-dimensional  wave  whereas  waves  with  a  wave  angle 
unequal  to  zero  are  strongly  amplified.  Since  the  resonance  mechanisms  of  a 
fundamental  or  a  subharmonic  resonance  is  dependent  on  the  amplitude  of 
the  primary  wave  and  not  on  the  amplitude  of  the  secondary  wave,  an  APG 
does  not  enhance  the  ’’classical”  fundamental  or  subharmonic  breakdown 
scenarios. 

Figure  6.24b  shows  the  lambda  leg  formation  typically  associated  with  this 
breakdown.  The  stretching  of  the  peak  stations  can  be  observed  in  the  span- 
wise  vorticity  isosurface  plot  (see  figure  6.24c).  Smaller  amplitudes  (<  1%  of 
the  free-stream  velocity)  of  the  primary  (1, 0)-wave  show  no  fundamental  res¬ 
onance  within  the  computational  domain  so  that  then  an  oblique  breakdown 
occurs. 

6.3.2  Classical  Subharmonic  Resonance 

While  Thumm  (1991)  was  able  to  find  a  weak  subharmonic  resonance  for 

Ma  =  1.6,  Bestek  &  Eissler  (1996)  could  not  successfully  find  it  for  Ma  =  4.8. 

At  Ma  =  3,  there  has  been  no  indication  of  a  subharmonic  resonance  with  a 

fundamental  2D  wave  in  our  simulations. 

* 

6.3.3  Oblique  Fundamental  Resonance 

To  find  a  breakdown  scenario  which  might  be  at  least  as  viable  as  the  oblique 
breakdown,  we  realized  that  a  fundamental  resonance  between  two  three- 
dimensional  waves  is  possible.  A  first  hint  was  provided  in  the  “classical”  fun¬ 
damental  (Klebanoff-type)  breakdown  where  the  (l,±2)-waves  are  strongly 
amplified  (c.f.  figure  6.23). 
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Figure  6.25  Oblique  fundamental  resonance:  Amplitude  distribution  of  u- 
pertubation  versus  downstream  direction.  Re  =  10s,  F  =  5  x  10-5,  = 

60°,*1>2  =  74°. 

Therefore,  (1,  ±1)  primary  disturbance  waves  (amplitudes  <  0.06%  of  the 
free-stream  velocity)  and  (1,  ±2)  secondary  disturbance  waves  (amplitudes  <. 
0.005%  of  the  free-stream  velocity)  are  combined  for  this  breakdown  scenario. 
Since  the  frequency  (F  =  5  x  10-5)  is  the  same  for  both  the  primary  and 
secondary  waves,  we  call  this  an  “oblique  fundamental  resonance”.  The  wave 
angles  were  chosen  to  be  60°  and  74°  for  the  primary  and  secondary  waves, 
respectively.  So  far,  for  the  parameters  tested,  we  found  that  the  oblique 
fundamental  resonance  in  comparison  with  the  “classical”  oblique  breakdown 
is  not  an  as  strong  mechanism. 

Although  figure  6.25  shows  a  deviation  from  the  linear  behavior  of  the  sec¬ 
ondary  wave,  the  amplitude  stays  about  15  times  smaller  than  the  primary 
wave.  Also  the  nonlinearly  generated  steady  (0,  l)-vortex  rises  to  the  same 
amplitude  level  as  the  secondary  wave,  while  the  (0,2)-  and  (l,3)-modes, 
associated  with  a  “classical”  oblique  breakdown  amplify  to  the  order  of  the 
primary  wave.  Because  the  flow  structures  shown  in  figures  6.26a  and  6.26b 
resemble  the  ones  of  an  oblique  breakdown,  we  believe  that  in  spite  of  the 
resonance  of  the  secondary  wave,  the  “classical”  oblique  breakdown  is  the 
superior  mechanism. 

6.3.4  Oblique  Subharmonic  Resonance 

Three  different  combinations  of  an  oblique  primary  disturbance  with  an 
oblique  secondary  disturbance  with  half  the  frequency  of  the  primary  wave 
are  possible.  Kosinov  &  Tumin  (1996)  concluded  that  a  subharmonic  reso¬ 
nance  of  two  three-dimensional  waves  might  play  an  important  role  in  super- 


Figure  6.26  Isosurface  of  Q-criterion  Q  =  0.2  (left)  and  spanwise  vorticity 
wz  =  15  (right)  for  a  oblique  fundamental  resonance;  Re  =  105,  Rx  =  705  — 
925,  F  =  5  x  IQ”5,  =  60°. 


Primary  wave 

Secondary  wave 

frequency 

wave  angle 

frequency 

wave  angle 

case  I 

F2i2  =  5  x  10"5 

$2,2  =  60° 

jFu  =  2.5  x  10-5 

=  60° 

case  II 

F2i1  =  5  x  10~5 

tt2,i  =  60° 

F1)2  =  2.5  x  10~6 

*1,2  =  82° 

case  III 

F2>  1  =  5  x  10"5 

^2,i  =  60° 

Fitl  =  2.5  x  IQ"5 

O 

II 

H 

Table  6.5  Oblique  Subharmonic  Breakdown.  Overview  of  studied  cases. 

sonic  boundary  transition  based  on  their  experimental  findings  at  Ma  =  2. 
Their  numerical  and  theoretical  work  confirmed  a  possible  resonance  between 
a  primary  (2,  l)-wave  with  secondary  (1,  -2)-  and  (1,  l)-waves.  In  addition, 
primary  (2, 2)-  and  secondary  (1,  l)-waves  can  both  travel  with  the  same 
wave  angle.  To  distinguish  these  subharmonic  resonances  from  the  “classi¬ 
cal”  one  which  involves  a  two  dimensional  primary  and  oblique  secondary 
disturbance,  we  call  this  breakdown  “oblique  subharmonic  resonance” . 

For  our  investigations  for  a  Mach  3  boundary  layer,  three  different  cases  were 
computed  (see  Table  6.5).  First,  the  primary  disturbances  (2,  ±2)  are  excited 
with  an  amplitude  of  0.3%  (of  the  free-stream  velocity),  and  the  secondary 
disturbances  (1,±1)  have  an  amplitude  of  0.005%.  Thus  both  primary  and 
secondary  waves  have  the  same  wave  angle  of  ^  =  60°.  For  case  II  and  case 
III,  secondary  waves  of  different  wave  angles  are  combined  with  a  primary 
wave  (2,  ±1)  (^2,1  —  60°,  amplitude  of  0.1%).  For  case  II,  the  secondary 
waves  are  (1,  ±2)-modes  with  $1>2  =  82°  and  also  disturbed  with  an  ampli¬ 
tude  of  0.005%.  For  case  III,  the  secondary  disturbances  are  (l,±l)-modes 
with  an  amplitude  of  0.005%  and  =  74°.  The  primary  disturbance  am¬ 
plitude  in  case  I  is  three  times  higher  than  in  the  other  two  cases  to  see  at 
least  a  small  deviation  from  the  linear  eigen-behavior  of  the  secondary  wave. 


Figure  6.27  Oblique  subharmonic  resonance:  Amplitude  distribution  of  u- 
perturbation  versus  downstream  direction  for  cases  I,  II  and  III  from  left  to 
right,  respectively;  Re  =  105,  F2>x  —  5  x  10-5,  FltX  =  2.5  x  10-5,  ^2,x  =  60°. 


Figure  6.28  Oblique  subharmonic  resonance:  Isosurface  of  Q-criterion  for 
cases  I  (Q  =  3),  II  (Q  =  0.2)  and  III  (Q  =  0.2)  from  left  to  right,  respectively; 
Re  =  105,  F2jX  =  5  x  1CT5,  Fhx  =  2.5  x  10"s,  tf2,x  =  60° 

But,  as  can  be  seen  in  figure  6.27a,  the  nonlinear  generation  of  higher  modes 
is  still  not  emphasized.  However,  the  amplification  curves  of  figures  6.27b 
and  6.27c  clearly  indicate  that  strong  nonlinear  interaction  can  occur  indi¬ 
cating  that  oblique  subharmonic  breakdown  mechanisms  may  be  relevant  for 
the  Mach  3  boundary  layer.  Surprisingly,  the  oblique  subharmonic  break¬ 
down  produces  a  more  rapid  amplitude  growth  in  the  downstream  direction 
at  lower  amplitudes  if  the  primary  and  the  secondary  waves  are  forced  with 
different  wave  angles  than  if  primary  and  secondary  waves  have  the  same 
wave  angle  (c.f.  figure  6.27a  with  figures  6.27b  and  6.27c).  This  is  an  in¬ 
dication  that  the  wave  angle  between  primary  and  secondary  waves  plays 
an  important  role  for  energy  transfer  from  the  base  flow  to  the  disturbance 
waves  and  between  primary  and  secondary  waves.  In  case  II,  where  the  sec¬ 
ondary  wave  is  mode  (1, 2),  the  nonlinearly  generated  (1,  l)-wave  also  shows 
some  nonlinear  resonance  slightly  downstream  of  the  resonance  location  of 
the  (1, 2)-wave.  But  in  case  III,  where  the  forced  (1,  l)-wave  shows  resonance 
about  the  same  location  as  in  case  II,  the  generated  (1, 2)-wave  reveals  no 
nonlinear  resonance. 

The  structures  in  figures  6.28a,  6.28b  and  6.28c  all  strongly  resemble  the 
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Figure  6.29  Oblique  subharmonic  resonance:  Isosurface  of  spanwise  vorticity 
uz  =  15  for  cases  I,  II  and  III  from  left  to  right,  respectively;  Re  =  10s, 
F2  x  =  5  x  1CT5,  FltX  =  2.5  x  10~5,  $2>a;  =  60° 

structures  found  in  an  oblique  breakdown  (see  figure  6.40  in  §  6.4).  For  the 
simulations  with  different  wave  angles  (case  II  and  case  III),  the  spanwise 
vorticity  (uz)  exhibits  an  additional  “tongue”  in-between  the  two  main  struc¬ 
tures  (see  figures  6.29b  and  6.29c  in  comparison  with  figure  6.26c).  Overall, 
it  appears  that  similar  mechanisms  as  for  the  “classical”  oblique  breakdown 
are  present. 

6.4  Oblique  Breakdown 

Thumm  (1991)  was  the  first  who  realized  the  importance  of  the  oblique 
breakdown  as  a  mechanism  that  could  lead  to  turbulence.  He  described  in 
detail  the  physics  of  the  start  up  of  the  oblique  breakdown.  Until  now  it  is 
still  not  clear  if  an  oblique  breakdown  really  leads  to  turbulence  and  if  this 
mechanism  plays  a  significant  role  in  nature. 


Figure  6.30  Wave  front  of  the  3D  waves  in  the  computational  domain. 


What  is  oblique  breakdown?  According  to  Mack  (1984),  the  Squire  theorem, 
which  states  that  the  critical  Reynolds  number  is  defined  by  the  location 
where  2D  waves  are  amplified,  is  not  valid  for  compressible  flow.  Instead,  3D 
waves  are  highly  amplified  and  for  a  given  Mach  number  there  is  a  specific 
wave  angle  which  leads  to  most  amplified  TS-waves.  Oblique  breakdown 
is  caused  by  the  nonlinear  interaction  of  these  3D  waves.  In  the  numerical 
simulations  presented  in  this  work,  the  development  of  a  3D  wave  throughout 
the  computational  domain  is  simulated,  similar  to  the  §  6.2,  but  with  higher 
amplitudes.  Since  the  DNS  code  employs  symmetric  Fourier  transforms  in 
the  spanwise  direction,  two  3D  waves,  with  #  and  — \F  as  their  wave  angles, 
are  excited  simultaneously  in  the  disturbance  slot  (figure  6.30).  If  the  3D 
waves  reach  large  amplitudes  farther  downstream,  they  will  start  to  interact 
with  themselves  and  with  each  other.  For  a  more  detailed  investigation  the 
flow  properties  are  decomposed  using  a  Fourier-transform  with  respect  to 
time  and  the  spanwise  direction: 

H  K 

<f>(x,y,z,t)  =  4>(x,y)ei(h^i+kyz\  (6.2) 

h=—H  k=—K 

where  H  is  the  number  of  Fourier  modes  in  time  and  K  in  the  ^-direction. 
This  formulation  has  the  advantage  that  it  can  be  used  to  show  how  higher 
modes  in  time  and  space  are  created  by  the  nonlinear  terms  of  the  Navier- 
Stokes  equations.  For  example, 


dui 


(6.3) 


creates  terms  like 


gi(hi0t+kiyz)  _  gilkipt+kiyz)  __  gi((hi+/i2)/3t+(fci +<12)72) 


(6.4) 


Equation  (6.4)  indicates  that  the  modes  in  time  and  in  the  2- direction  of 
the  interacting  waves  must  be  added  to  generate  higher  solutions.  As  an 
abbreviation  equation  (6.4)  can  also  be  written  as 


(fill  &i)  +  (fi.2)  ^2)  =  (fix  +  fi2,  fii  +  £2);  (6.5) 

where  (fii,&i)  indicates  a  mode  combination  with  the  frequency  hi/3  and 
the  spanwise  wave  number  A77.  For  a  3D  wave  with  a  frequency  j3  and  a 
wave  number  7,  fi  and  k  are  equal  to  one.  The  wave,  which  is  mirrored  at  the 
symmetry  line  of  the  domain,  which  is  parallel  to  x  and  located  in  the  middle 
of  the  spanwise  direction,  has  — as  its  wave  angle  and,  consequently,  fi  =  1 
and  k  =  —1.  The  possible  mode  combinations,  that  are  created  by  these  two 
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waves,  are  listed  in  table  6.6.  The  resulting  waves  exist  in  addition  to  the 
original  waves  and  can  also  interact  with  them  to  produce  even  higher  mode 
combinations  like  (1,3)  or  (4,0).  An  interesting  consequence  of  the  oblique 
breakdown  is  that  mode  combinations  with  odd  modes  in  time  have  only- 
odd  modes  in  z  and  with  even  modes  in  time  only  even  modes  in  2  (Thumm, 
1991). 


(/ll,  &l)  +  (/l2,  £2) 

(hi  +  /1-2,  ki  +  £2) 

(1,1) +(1,1) 

(2,2) 

(1,1) +  (1,-1) 

(2,0) 

(1, 1)-+  (-1,1) 

(0,2) 

(1,1)  +  (-1,-1) 

(0,0) 

(1,  —1)  +  (1,  —1) 

(2,-2) 

(1,  -1)  +  (—1,  -1) 

(0,-2) 

(-1>  1)  +  (-1,  1) 

(-2,2) 

(-1,1) +  (-1,-1) 

(-2,0) 

(-1,-1)  +  (-1,-1) 

(-2,-2) 

Table  6.6  Selected  mode  combinations  of  a  (1,1)  and  (1,-1)  wave-pair. 


Figure  6.31  Wave  number  diagram  for  the  mode  combinations  (1,1),  (1,-1), 
(2,0)  and  (0,2). 

The  mode  combinations  can  be  plotted  as  vectors  in  a  wave  number  diagram 
(figure  6.31).  The  generation  of  the  higher  mode  combinations  is  equiva¬ 
lent  to  a  vector  summation  or  subtraction.  In  figure  6.31  two  examples  are 


shown.  Note  that  (1,1)  and  (1,-1)  are  the  subharmonic  mode  combinations 
of  (2,0).  However,  there  is  no  subharmonic  resonance  (staggered  A- vortices) 
since  the  3D  waves  are  linearly  more  amplified  than  the  2D  wave.  Thumm 
(1991)  observed  that  the  mode  combination  (0,2)  has  the  highest  influence 
on  oblique  breakdown  in  its  start  up. 

It  is  of  interest  to  see  how  a  mode  combination  would  look  in  physical  space. 
Since  we  are  interested  in  vorticity  waves,  i.e.,  TS-waves,  a  wave  ansatz  is 
introduced  which  has  the  form 

y,  Z,  t )  =  foyfrax+kyz+hact)'  (6.6) 

If  h  is  increased  the  wave  will  have  a  smaller  wavelength  (divided  by  h) 
and  if  k  is  increased  the  wave  will  have  a  larger  wave  angle.  <p  represents  a 
arbitrary  flow  quantity,  e.g.,  the  u- velocity.  An  instantaneous  plot  of  the  u- 
velocity  of  the  wave  at  a  specific  y- location  includes  positions  of  the  maximum 
and  minimum  of  the  u- velocity.  The  positions  of  the  maximum  (minimum) 
are  repeated  every  wavelength  as  is  illustrated  in  figure  6.32.  The  plot  is 
calculated  with  the  real  part  of  equation  (6.6),  but  without  amplification  of 
the  waves.  The  values  on  the  2-axis  do  not  match  any  values  obtained  from 
a  DNS,  since  they  are  intended  for  visualization  purposes  only.  Since  the 


Figure  6.32  Example  of  a  3D  wave  traveling  with  the  wave  angle  T  =  65°. 

u- velocity  vector  is  parallel  to  the  flow  direction,  a  wave  in  the  u- velocity  is 
a  longitudinal  wave. 

The  interference  of  waves  with  the  wave  angle  T  and  — $  is  visualized  in 
figure  6.33.  It  represents  a  superposition  of  two  3D  waves.  The  linear  stage 
of  oblique  breakdown  should  be  dominated  by  such  structures.  They  have  a 
shape  similar  to  the  letter  “X”  where  the  maxima  are  located  in  the  center 
and  at  the  ends  of  the  “X”.  Structures  with  a  “X” -shape  in  the  pressure 
distribution  were  observed  in  preliminary  simulations  of  the  nonlinear  stage 


of  oblique  breakdown,  and  we  will  look  for  them  in  the  results  presented  in 
the  following  sections.  Before  the  results  obtained  by  the  DNS  are  discussed, 
attention  is  being  devoted  to  the  influence  of  the  (0,2)  and  (0,4)  mode  combi¬ 
nation  on  the  “X” -shape.  The  (0,2)  and  the  (0,4)  are  stationary  distortions 
of  the  baseflow.  Figures  6.34  and  6.35  illustrate  how  the  distortion  of  the 
baseflow  alters  the  “X” -Shape. 


Figure  6.35  Influence  of  mode  combination  (0,4). 
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6.4.1  Without  Pressure  Gradient 

All  oblique  breakdown  calculations  are  conducted  with  the  flow  properties 
which  were  used  for  the  supersonic  validation  cases  in  §  4.5.  A  detailed 
summary  of  the  computational  parameters  for  case  OBLNP3,  the  simulation 
presented  in  this  section,  can  be  found  in  appendix  C.  The  computational 
domain  starts  at  Rx  —  750  and  ends  at  Rx  =  1438.  Disturbances  with  the 
amplitude  Ait\  =  0.003  and  angular  frequency  F  =  3  x  10~5  are  introduced 
into  the  computational  domain  between  Rx  =  927  and  Rx  =  993.  The 
computational  domain  is  four  boundary  layer  thicknesses  high.  The  buffer 
domain  starts  at  Rx  —  1380.  At  the  location  where  the  disturbances  are 
introduced,  the  grid  is  equidistant  in  x  and  22  points  resolve  one  streamwise 
wavelength  Xx.  Downstream  of  the  disturbance  slot,  the  grid  is  stretched  until 
the  position  Rx  =  1365.  At  this  location  Ax  has  such  a  value  that  176  points 
resolve  one  streamwise  wavelength  Xx.  55  modes  are  used  in  the  spanwise 
direction.  The  resolution  of  the  calculation  is  higher  than  the  resolution  of 


Figure  6.36  Isosurface  of  streamwise  vorticity  \ux\  =  \  between  Rx  =  1085 
and  Rx  =  1291. 

the  simulations  presented  by  Thumm  (1991)  and  Eissler  (1995).  For  that 
reason  it  is  possible  to  compute  deeper  into  transition  and  to  investigate 
how  new  structures  arise.  Thumm  (1991),  in  his  work,  discussed  structures 
visualized  by  the  streamwise  vorticity  ux  which  have  the  shape  of  tongues. 
The  tongues  split  into  two  tips  further  downstream.  Preliminary  simulations 
carried  out  for  the  present  research  had  shown  that  for  higher  Mach  numbers 
the  tips  developed  more  strongly  on  one  side.  In  figure  6.36  these  structures 
are  plotted  for  \ux\  =  1  between  Rx  =  1085  and  Rx  =  1291  for  the  simulation 
OBLNP3.  The  top  view  of  figure  6.36  also  shows  the  previously  discussed 
“X” -shape.  However,  the  structures  do  not  really  look  like  the  predicted 
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Figure  6.37  Contour  of  u-velocity  between  Rx  =  1085  and  Rx  =  1291  at 
location  y  =  0.001  and  y  =  003. 


shape  shown  in  figure  6.33.  For  that  reason,  the  u-velocity  was  checked 
whether  there  the  “X-shape”  appears.  In  figure  6.37  contours  of  the  u- 
velocity  at  different  y-locations  in  the  boundary  layer  are  plotted.  At  y  = 
0.003  the  predicted  “X” -shape  is  visible.  Upstream,  two  parallel  black  lines 
are  getting  stronger  which  are  caused  by  the  (0,2)  mode  combination.  It 


Figure  6.38  Distribution  of  |uIOj]U*  (left)  and  u[1>k]max  (right)  versus  Rx. 

is  very  interesting  that  at  a  different  y-location  the  contours  look  differently. 
At  y  —  0.001  the  “X”  is  almost  not  visible  anymore  and  only  two  parallel 
black  lines  (the  (0,2)  mode  combination)  dominate  the  velocity  field.  The 
explanation  for  that  behavior  can  be  found  in  the  following  two  figures.  The 
mode  combinations  (0,2)  and  (1,1)  are  strongly  amplified  and  they  have  high 
amplitudes  between  Rx  =  1085  and  Rx  =  1291  (figure  6.38).  The  amplitude 
distribution  of  these  two  mode  combinations  over  y  for  different  downstream 


Figure  6.39  Development  of  the  amplitude  distribution  of  U[h,k}/'U[hji]ma «  over 
V- 

locations  is  illustrated  in  figure  6.39.  The  maximum  of  the  (1,1)  is  located 
at  a  different  location  than  the  maximum  of  the  (0,2).  This  means  that  the 
mode  combinations  have  a  different  influence  for  changing  y  locations.  At 
y  =  0.001  the  (1,1)  is  very  weak  and  that  is  why  the  (0,2)  is  strongly  visible 
in  figure  6.37.  At  y  =  0.003  both  modes  are  present.  It  is  also  of  interest 
that  the  (0,2)  mode  combination  is  significantly  changing  in  the  downstream 
direction.  Its  maximum  is  shifted  to  higher  values  of  y  and  the  absolute  value 
of  the  minimum  is  decreasing  until  position  R*  =  1357  where  the  locations 
of  the  maximum  and  the  minimum  are  inverted.  Note  that  at  Rx  =  1085 
the  absolute  value  of  the  minimum  is  higher  than  the  maximum.  At  position 
Rx  =  1208  the  maximum  is  higher  than  the  absolute  value  of  the  minimum. 
For  that  reason  there  is  a  significant  change  in  the  slope  of  the  amplitude 
distribution  of  the  (0,2)  at  around  Rx  —  1140  in  figure  6.38.  The  shift  of 
the  maximum  to  higher  y- values  is  the  explanation  for  the  changing  pattern 
of  the  u-velocity  contour  at  location  y  =  0.003  (figure  6.37).  Upstream,  the 
flow  is  dominated  by  the  “X”-shape.  Farther  downstream,  the  (0,2)  is  getting 
stronger  because  its  amplitude  value  is  changing  from  a  negative  to  a  positive 
value  for  higher  local  Reynolds  numbers  (until  Rx  =  1357). 

The  latter  discussion  emphasizes  that  in  the  u- velocity,  the  predicted  “X”- 
shape  can  be  found  since  the  “X” -shape  is  caused  by  the  mode  combinations 
(1,1)  and  (1,-1).  However,  it  is  not  clear  if  this  is  a  dominant  structure  since 
in  every  quantity  different  patterns  arise  at  different  y-locations.  A  good 
identification  of  structures  can  be  obtained  by  employing  a  vortex  identifica¬ 
tion  technique  called  Q- criterion.  It  is  based  on  the  vortex  definition  given 
by  Chong  et  al.  (1990).  A  vortex  is  located  in  a  region  in  the  flow  field 
where  the  vorticity  is  sufficiently  strong  to  cause  the  velocity  gradient  tensor 
to  be  dominated  by  the  rotation  tensor.  That  means  mathematically  that 
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Figure  6.40  “X”-shape  vortices  identified  by  Q-criterion  with  Q  =  10  between 
Rx  =  1085  and  Rx  =  1291. 


Figure  6.41  “X”-shape  vortices  identified  by  Q-criterion  with  Q  =  30  between 
Rx  =  1085  and  Rx  =  1291. 

the  magnitude  of  the  symmetric  part  of  the  velocity  gradient  tensor,  which 
describes  the  strain,  is  smaller  than  the  asymmetric  part  that  describes  the 
rotation.  Q  identifies  the  difference  between  the  magnitude  of  rotation  and 
the  magnitude  of  strain,  so  that 

Q  =  \(WiiWij-SijSij)i  (6.7) 

where  W{j  indicates  an  element  of  the  rotation  tensor  and  Sq  an  element 
of  the  strain  tensor.  For  positive  values  of  <5,  rotation  is  stronger  than 
strain.  Note  that  equation  (6.7)  is  an  incompressible  definition.  For  an 
incompressible  flow,  Q  is  invariant  with  respect  to  its  reference  frame.  That 


is  not  true  for  compressible  flows.  Nevertheless,  it  still  defines  locations 
where  rotation  dominates  the  flow  field.  If  equation  (6.7)  is  applied  to  the 
previously  discussed  flow  field  it  visualizes  structures  presented  in  the  figures 
6.40  and  6.41.  The  Q-criterion  confirms  the  “X”-shape. 

Since  the  integration  domain  extended  further  downstream  than  at  Rx  = 
1291,  the  deformation  of  the  “X” -shape  can- be  investigated.  Figures  6.38 
and  6.39  already  give  an  impression  of  what  will  happen  further  downstream. 
Higher  mode  combinations  are  gaining  importance.  The  amplitude  distribu¬ 
tion  in  y  of  the  (0,2)  (figure  6.39)  is  significantly  changing  between  Rx  =  1333 
and  Rx  =  1357.  This  is  also  connected  with  its  amplitude  distribution  in 
x  where  its  value  is  strongly  decreasing  downstream  of  Rx  =  1320  (figure 
6.39).  The  mode  combination  (0,4)  has  the  highest  value  of  all  stationary 
mode  combinations  downstream  of  Rx  =  1300.  Its  effect  on  the  w- velocity 
is  demonstrated  in  figure  6.42.  The  four  bright  lines  parallel  to  the  x-axis 
represent  the  four  maxima  of  the  (0,4).  At  the  end  of  the  domain  six  black 
lines  are  getting  stronger  which  are  the  minima  of  the  (0,6).  A  3D  view  of 


Figure  6.42  u- velocity  at  y  =  0.001  between  Rx  =  1333  and  Rx  =  1380. 

an  isosurface  of  the  u- velocity  for  the  value  u  —  0.65  is  given  in  figure  6.43. 
The  previously  discussed  “X” -shape  is  visible.  Further  downstream,  the  legs 
of  the  “X”  which  are  closer  to  the  buffer  domain  are  pushed  away  from  the 
wall.  This  process  clearly  thickens  the  boundary  layer  (figure  6.44).  Parti¬ 
cles  with  high  velocity  are  transported  closer  to  the  wall  and  particles  with 
low  velocity  are  pushed  away  from  the  wall.  At  the  end  of  the  domain  in 
figure  6.43  new  structures  arise  parallel  to  the  x-axis  at  z=  \z/2  (along  the 
symmetry  line).  The  Q-criterion  confirms  the  new  structure  in  figure  6.45. 
At  the  start  of  the  buffer  domain  (Rx  =  1380),  the  resolution  becomes  too 
coarse  to  resolve  the  small  structures.  Grid  mesh  due  to  the  lack  of  resolution 


Figure  6.43  3D  isosurface  of  the  u- velocity  with  u  =  0.65  between  Rx  — 1208 
and  Rx  =  1380. 


iBiunm 


X 


mabgp  siSSn  KPtertio 


taW  -txton 

Figure  6.44  Contours  of  the  u- velocity  in  the  computational  domain  between 
Rx  =  1208  and  Rx  =  1380  at  z  =  0.273  x  \z,  the  position  where  one  leg  of 
the  “X”  is  pushed  away  from  the  wall. 


occurs  in  the  ^-direction  (figure  6.45).  The  influence  of  higher  modes  in  time 
on  the  generation  of  the  new  structures  is  limited.  Only  the  mode  combina¬ 
tions  (2,0),  (2,2),  (2,4)  and  (2,6)  have  high  enough  amplitudes  (figure  6.46) 
to  significantly  alter  the  flow  field. 

6.4.2  With  Adverse  Pressure  Gradient 

In  this  section  two  simulations  are  compared,  one  with  (case  OBLPG)  and 
one  without  pressure  gradient  (case  OBLNP4).  The  computational  domain 
and  flow  quantities  are  the  same  with  the  exception  of  the  pressure  gra¬ 
dient  and  the  disturbance  amplitude  which  is  almost  four  times  higher  for 
the  simulation  without  APG.  The  domain  starts  at  Rx  =  1226  and  ends  at 
Rx  =  1782.  Disturbances  with  the  angular  frequency  of  F  =  3  x  10~5  are  in¬ 
troduced  into  the  domain  through  a  blowing  and  suction  slot  located  between 


Figure  6.45  “New  structures”  between  the  legs  of  the  “X” -Shape  identified 
by  the  Q-criterion  for  Q  =  500  between  Rx  =  1333  and  Rx  =  1380. 


Figure  6.46  Distribution  of  up^max  (left)  and  u^max  (right)  versus  Rx. 


Rx  =  1332  and  Rx  =  1369.  The  buffer  domain  starts  at  Rx  =  1715.  The 
domain  is  about  four  boundary  layer  thicknesses  high.  The  grid  stretching 
starts  downstream  of  the  disturbance  slot  and  ends  at  Rx  =  1657.  The  ratio 
is  chosen  in  such  way  that  88  points  are  in  one  streamwise  wavelength  Xx 
downstream  of  Rx  =  1657.  16  modes  are  used  to  resolve  the  spanwise  direc¬ 
tion.  With  this  low  resolution  only  the  start  up  of  oblique  breakdown  with 
and  without  APG  can  be  simulated.  It  was  not  possible  to  compute  farther 
into  transition  and,  therefore,  to  choose  a  higher  resolution  since  the  pressure 
distribution  at  the  free-stream  drifts  from  its  desired  value.  This  may  well  be 
the  correct  physical  behavior;  however,  more  investigations  are  required  to 
prove  this.  The  simulation  OBLPG  has  a  moderate  APG.  The  distribution 
of  the  Hartree  parameter  (3h  over  Rx  indicates  an  almost  constant  value  of 
/ 3h  «  —0.12  (figure  6.47). 
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Figure  6.47  Distribution  of  Hartree  parameter  f3H  versus  R*  for  case 
OBLNP4. 


Figure  6.48  Distribution  of  |u[0,k]U*  (top)  and  u[1>k]max  (bottom)  versus  R* 
for  cases  OBLNP4  (no  APG,  left)  and  OBLPG  (with  APG,  right). 

The  amplitude  distributions  of  the  individual  mode  combinations  (h,k)  versus 
Rx  for  both  simulations  are  plotted  in  figure  6.48.  All  mode  combinations 
(h,k)  experience  a  much  higher  amplification  for  case  OBLPG  than  for  case 
OBLNP4.  The  APG  has  a  strong  impact  on  the  amplification  rate.  For  that 
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reason  the  disturbance  amplitude  in  the  suction  and  blowing  slot  has  to  be 
much  smaller  than  for  case  OBLNP4,  as  mentioned  earlier,  otherwise  the 
disturbance  waves  would  reach  much  faster  the  nonlinear  stage  and  the  flow 
would  start  to  transition  further  upstream.  However,  this  is  not  desired  since 
in  both  simulations  the  flow  is  required  to  transition  shortly  upstream  of  the 
buffer  domain,  so  that  the  calculation  can  run  stable  and  can  converge  out. 
The  abrupt  changes  of  the  slope  of  the  stationary  mode  combinations  (0,2) 
and  (0,4)  in  figure  6.48  (top  left)  can  be  explained  in  the  same  way  like  for 
figure  6.38.  For  example,  the  absolute  value  of  the  minimum  (of  the  (0,2))  is 
at  one  JV  location  bigger  and  at  another  Rx  smaller  than  the  absolute  value 
of  the  maximum  (also  of  the  (0,2)),  so  that  at  the  position  where  it  changes  an 
abrupt  switch  in  the  slope  of  the  amplitude  distribution  over  Rx  is  visible. 
In  the  case  OBLPG  abrupt  changes  cannot  be  found  in  figure  6.48  (top 
right).  However,  they  exist  at  lower  amplitude  levels  and  therefore  at  lower 
local  Reynolds  numbers.  The  mode  combination  (1,1)  for  case  OBLNP4  is 
weakly  amplified  (figure  6.48,  bottom  left).  The  amplitude  value  of  the  (1,3) 
is  increasing  much  faster  than  the  amplitude  value  of  the  (1,1).  For  that 
reason,  downstream  of  Rx  «  1650  the  (1,3)  has  a  higher  amplitude  than  the 
(1,1).  For  the  case  OBLPG,  the  (1,1)  is  strongly  amplified  and  that  is  the 
reason  why  the  (1,3)  does  not  reach  higher  amplitude  values  than  the  (1,1) 
in  the  presented  computational  domain. 


Figure  6.49  Isosurface  of  Q  =  2  for  case  OBLNP4  between  Rx  =  1644  and 
Rx  =  1715. 


The  structures  which  are  identified  by  the  Q-criterion  look  very  similar  to 
the  structures  obtained  by  simulation  OBLNP3  (figure  6.49,  6.50,  6.51  and 
6.52).  The  “X” -shape  is  dominant  for  both  calculations.  In  the  simulation 
OBLNP4,  the  maximum  in  the  center  of  the  “X”  is  divided  in  three  parts  and 
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Figure  6.50  Isosurface  of  Q  =  5  for  case  OBLPG  between  Rx  =  1644  and 
Rx  =  1715. 


Figure  6.51  Isosurface  of  Q  =  30  for  case  OBLNP4  between  Rx  =  1644  and 
Rx  =  1715. 


two  spikes  are  coming  out  of  its  head  for  small  values  of  Q.  This  behavior  was 
also  observed  in  simulation  OBLNP3  (however,  no  plot  in  the  last  chapter 
indicated  this  because  attention  was  being  devoted  to  different  properties 
of  the  “X” -shape).  Most  likely  the  division  of  the  maximum  is  caused  by 
the  (1,3)  mode  combination  since  the  calculation  with  APG  does  not  show 
this.  Nevertheless,  detailed  investigations  have  to  be  performed  to  prove  this 
statement.  The  slim  structure  along  the  symmetry  line  in  figure  6.50  at  the 
end  of  the  domain  can  also  be  observed  in  the  simulation  without  APG  at  a 
lower  level  of  Q.  It  does  not  serve  as  an  indication  for  the  differences  between 
both  calculations. 
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Figure  6.52  Isosurface  of  Q  =  50  for  case  OBLPG  between  Rx  =  1644  and 
Rx  =  1715. 


6.5  Summary  of  Nonlinear  Breakdown  Scenarios 

An  overview  of  the  nonlinear  breakdown  scenarios  we  investigated  can  be 
found  in  table  6.7.  As  already  discussed  in  §  2.3,  in  addition  to  the  clas¬ 
sical  subharmonic  and  fundamental  breakdown  scenarios  of  incompressible 
boundary  layers,  the  so-called  “oblique  breakdown”  may  play  a  role  for  su¬ 
personic  boundary  layers.  For  the  oblique  breakdown,  only  modes  (1,±1) 
have  to  be  forced.  We  first  discovered  the  oblique  breakdown  mechanism 
for  Ma  =  1.6  (Thumm,  1991),  and  it  was  found  later  that  it  may  also  be 
relevant  for  Ma  =  4.8  (Eissler,  1995).  Thus,  we  first  needed  to  confirm  if 
the  oblique  breakdown  was  also  relevant  for  the  Ma  =  3  case  of  the  Prince¬ 
ton  experiments.  Nonlinear  transition  scenarios  for  Ma  =  3  are  likely  to  be 
more  interesting  than  for  Ma  =  1.6  (albeit  also  more  complicated)  because, 
depending  on  the  local  Reynolds  number,  first  and  second  modes  can  both 
be  excited  linearly  and  nonlinearly. 

In  flight  vehicles,  boundary  layers  are  subject  to  streamwise  pressure  gradi¬ 
ents,  both  favorable  and  adverse.  Adverse  pressure  gradients  are  more  critical 
as  they  accelerate  transition  while  favorable  pressure  gradients  delay  tran¬ 
sition.  Therefore,  in  anticipation  of  the  planned  experiments  with  adverse 
pressure  gradients  at  the  Princeton  facility,  we  started  oblique  breakdown 
simulations  (under  the  conditions  of  the  Princeton  experiments)  to  investi¬ 
gate  if  the  oblique  breakdown  is  also  relevant  for  Ma  =  3  when  an  adverse 
pressure  gradient  is  imposed.  As  is  characteristic  for  the  oblique  breakdown, 
only  (1,  ±1)  modes  are  forced,  in  this  case  two  oblique  waves  with  a  wave  an¬ 
gle  of  65  degrees  (most  amplified  according  to  linear  theory).  For  the  forcing, 
a  blowing  and  suction  slot  is  used  that  is  located  one  streamwise  wavelength 
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Table  6.7  Overview  of  breakdown  scenarios:  (h,  ±k,  where  h  is  the  frequency, 
k  the  spanwise  Fourier  mode;  ■  primary  wave,  Q  secondary  wave. 


(A*)  downstream  of  the  inflow  boundary.  It  was  shown  that  for  the  adverse 
pressure  gradient  case  (Hartree  parameter  fix  varied  between  -0.11  and  -0.12 
in  the  computational  domain)  the  nonlinear  wave  components  grow  more 
rapidly  than  for  the  zero  pressure  gradient  case,  a  manifestation  that  an  ad¬ 
verse  pressure  gradient  does  not  only  not  interfere  with  this  mechanism  but 
rather  enhances  it. 

To  study  the  mechanisms  present  in  an  oblique  breakdown  in  more  detail, 
a  high  resolution  simulation  which  extended  farther  into  the  late  stages  of 
transition  was  also  investigated.  The  oblique  modes  (1,  ±1)  were  rapidly 
amplified  in  the  downstream  direction  and,  as  their  amplitudes  grew,  mul¬ 
tiple  other  modes  were  generated  from  nonlinear  interactions.  Prominent  in 
this  respect  was  the  generation  of  steady  longitudinal  modes  [(0,2),  (0,4), 
...],  which  are  steady  vortices.  This  implies,  that  steady  vortices  (such  as, 
for  example,  those  produced  by  distributed  roughness)  could  also  trigger  this 
breakdown  mechanism  when  only  very  small  oblique  disturbances  are  present. 
Typical  flow  structures  that  arise  during  the  later  stages  of  this  oblique  break¬ 
down  scenario  were  documented  in  §  6.4.  An  important  observation  gained 
from  these  simulations  is  the  fact  that  the  nonlinear  mechanism  of  the  oblique 
breakdown  sets  in  at  very  low  amplitudes  (=  0.3%  of  the  free-stream  velocity) 
and  thus  is  a  potentially  highly  relevant  breakdown  mechanism  for  Ma  —  3 
at  low  free-stream  turbulence  levels  (free  flight).  We  therefore  recommend 
strongly  future  investigations  of  this  mechanism,  preferably  in  collaboration 
with  the  Princeton  experiments. 

In  addition  to  the  oblique  breakdown,  we  investigated,  if  the  “classical”  fun¬ 
damental  (Klebanoff)  breakdown  scenario  [(1,0),  (1,  ±1)]  is  also  relevant. 
The  amplification  curves  and  the  structures  resulting  from  the  “classical” 
fundamental  breakdown  were  shown  in  §  6.3.  However,  since  much  larger 
amplitudes  were  required  to  initiate  the  fundamental  breakdown,  it  is  likely 
to  be  relevant  only  for  high  free-stream  turbulence  conditions.  Nevertheless, 
this  breakdown  mechanism  needs  to  be  investigated  because  of  its  possi¬ 
ble  relevance  for  interpreting  experimental  measurements  from  “noisy”  (high 
free-stream  turbulence  level)  facilities.  It  should  be  noted  that  measurements 
in  the  Russian  experiments  for  a  sharp- nosed  cone  at  Ma  =  6  (Shiplyuk  et  ai, 
2003)  indicated  that  a  two-dimensional  second-mode  subharmonic  instabil¬ 
ity  was  present.  In  our  Navier-Stokes  simulations  at  Ma  =  3,  resonance  oc¬ 
curred  at  much  lower  amplitudes  for  the  second-mode  “classical”  fundamental 
breakdown  scenario  than  for  the  first-mode  scenario.  Therefore,  “classical” 
fundamental  and  “classical”  subharmonic  breakdown  scenarios  are  likely  im¬ 
portant  for  second-mode  disturbances. 

In  addition  to  the  “classical”  nonlinear  resonance  mechanisms,  we  investi¬ 
gated  other  resonance  mechanisms  that  may  become  relevant  due  to  the  fact 


that,  for  supersonic  boundary  layers,  the  first-mode  oblique  disturbances  are 
linearly  more  amplified  than  first-mode  two-dimensional  disturbances  (which 
are  most  amplified  in  incompressible  flows).  Assuming  that  this  behavior  also 
holds  for  the  nonlinear  stages,  we  investigated  subharmonic  and  fundamen¬ 
tal  nonlinear  resonance  mechanisms  for  oblique  fundamental  waves  where 
both  primary  and  secondary  disturbance  waves  are  three-dimensional.  For 
the  oblique  fundamental  resonance,  disturbances  in  the  (1,±1)  modes  for  the 
primary  waves  (amplitudes  =  0.02%  of  the  free-stream  velocity)  and  in  the 
(1,±2)  modes  for  the  secondary  waves  (amplitudes  =  0.005%  of  the  free- 
stream  velocity)  resulted  in  a  wave  angle  of  60  degrees  and  74  degrees  for 
the  primary  and  secondary  waves,  respectively.  Our  results  indicate  that 
the  oblique  fundamental  breakdown  triggers  an  essentially  “oblique  break¬ 
down,”  however,  at  even  lower  disturbance  amplitudes  than  when  only  one 
pair  of  oblique  waves  is  introduced.  Simulations  with  different  free-stream 
temperatures  showed  that  the  influence  of  the  free-stream  temperature  on 
the  nonlinear  development  is  not  profound  as  the  onset  of  transition  moves 
only  slightly  upstream  for  T  =  103.6^  when  compared  with  T  =  300A.  The 
extent  of  the  nonlinear  transition  regime  is  somewhat  larger  for  T  =  103.6AT 
than  for  T  =  300AT.  Therefore,  an  outflow  treatment  (buffer  domain)  had  to 
be  used  for  this  calculation  in  order  to  prevent  a  “blow  up”  in  the  case  with 
T  =  300 K  with  equal  resolution. 

In  §  6.3,  we  also  investigated  the  subharmonic  resonance  that  Kosinov  &  Tu- 
min  (1996)  suggested  to  be  relevant  in  the  experiments  of  a  Mach  2  bound¬ 
ary  layer.  Three  different  scenarios  were  considered.  First,  the  primary 
disturbances  (2, ±2)  were  excited  with  an  amplitude  of  0.3%  (of  the  free- 
stream  velocity),  and  the  secondary  disturbances  (1,±1)  had  an  amplitude 
of  0.005%.  Thus  both  primary  and  secondary  waves  had  the  same  wave 
angle  of  ^  =  60  degrees.  For  the  second  and  the  third  case,  secondary 
waves  of  different  wave  angles  were  combined  with  a  primary  wave  (2,±1) 
(i>  =  60  degrees,  amplitude  of  0.1%).  For  the  second  case,  the  secondary 
waves  were  (1,±2)  modes  with  ip  =  81  degrees  and  also  disturbed  with  an 
amplitude  of  0.005%.  For  the  third  case,  the  secondary  disturbances  were 
(1,±1)  modes  with  an  amplitude  of  0.005%  and  ip  =  74  degrees.  The  ampli¬ 
fication  curves  clearly  indicated  that  strong  nonlinear  interaction  can  occur. 
Surprisingly,  the  oblique  subharmonic  breakdown  produced  a  more  rapid 
amplitude  growth  in  the  downstream  direction  at  lower  amplitudes  if  the 
primary  and  the  secondary  waves  are  forced  with  different  wave  angles  than 
if  primary  and  secondary  waves  travel  with  the  same  wave  angle.  This  is  an 
indication  that  the  wave  angle  between  primary  and  secondary  waves  plays 
an  important  role  for  energy  transfer  from  the  baseflow  to  the  disturbance 
waves.  The  different  wave  angles  between  primary  and  secondary  waves  also 


have  an  effect  on  the  flow  structures  during  breakdown,  in  particular  also  on 
the  spanwise  vorticity  structures.  Again,  the  experiments  in  the  Princeton 
tunnel  and  more  detailed  numerical  investigations  need  to  confirm  this. 


I)  Oblique  bmdulown 


II)  ‘Cki&sJcur  fundamental  breukdown 


Figure  6.53  Near-wall  temperature  distribution  for  different  breakdown  sce¬ 
narios.  I)  Oblique  breakdown,  F  =  3  x  10-5,  =  103 .6K.  a)  without 

pressure  gradient.  b)with  pressure  gradient.  II)  “Classical”  fundamental 
breakdown,  1st  mode,  F  =  5  x  10-5,  =  103.6AT.  Ill)  Oblique  fundamen¬ 

tal  breakdown,  F  =  2  x  10"5.  a)  TTO  =  103.6A.  b)  =  300AT.  IV)  Oblique 
subharmonic  breakdown,  F2> i  =  5  x  10-5,  Fi>2  =  2.5  x  10-5,  =  103.6AT. 
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As  can  be  seen  in  figure  6.53,  the  various  breakdown  scenarios  also  leave  a 
strong  imprint  on  the  near- wall  temperature  distribution.  This  is  of  signifi¬ 
cance,  of  course,  as  heat  loads  associated  with  laminar-turbulent  transition 
are  a  major  concern  in  the  design  and  safe  operation  of  high-speed  flight  vehi¬ 
cles.  Interesting  in  this  respect  is  the  fact  that  certain  breakdown  scenarios 
such  as  the  oblique  and  oblique  subharmonic  resulted  in  higher  tempera¬ 
tures  than  the  various  fundamental  breakdown  scenarios  with  the  “classical” 
(Klebanoff)  breakdown  resulting  in  the  lowest  temperatures.  This  may  possi¬ 
bly  be  exploited  by  alternating  the  natural  transition  process  using  artificial 
means,  such  as  surface  texturing,  or  changing  the  design  altogether  so  that 
certain  disadvantageous  (from  a  wall  heating  point  of  view)  breakdown  sce¬ 
narios  are  avoided.  Prom  an  experimental  point  of  view,  it  is  of  interest  if 
the  various  nonlinear  breakdown  scenarios  could  be  observed  with  visual¬ 
ization  techniques  that  are  realizable  for  supersonic  flows,  such  as  Schlieren 
techniques.  Therefore,  in  figure  6.54,  Schlieren  pictures  obtained  from  the 
various  breakdown  simulations  are  shown.  Clearly,  the  different  breakdown 
mechanisms  leave  a  distinct  mark  on  the  Schlieren  patterns.  Notice,  for  ex¬ 
ample,  the  rope-like  structures  that  were  observed  in  the  experiments  (see 
Pruett  &  Chang,  1995).  Therefore,  when  simulations  are  combined  with 
experiments,  dominant  breakdown  mechanisms  may  be  identified  from  ex¬ 
perimental  Schlieren  pictures  before  going  into  detailed  and  expensive  mea¬ 
surements. 


side  view  top  view 

I  a)  II  a) 


Figure  6.54  Simulated  Schlieren  pictures  of  various  breakdown  mechanisms. 
Shown  are  side  view  of  wall-normal  density  gradient  (I)  and  top  view  of 
streamwise  density  gradient  (II):  a)  oblique  breakdown  with  zero  pressure- 
gradient,  b)  oblique  breakdown  with  APG,  c)  classical  fundamental  break¬ 
down,  d)  oblique  fundamental  breakdown,  e)  oblique  subharmonic  break¬ 
down. 
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7.  CONCLUSIONS 

With  the  funding  from  the  current  AFSOR  grant,  we  focused  our  numer¬ 
ical  simulation  efforts  on  a  collaboration  with  the  experiments  for  a  Mach 
3  flat-plate  boundary  layer  carried  out  at  Princeton  University  (G.  Brown). 
Brown  and  co-workers  (see  Graziosi,  1999;  Graziosi  k  Brown,  2002)  have 
been  performing  “natural”  and  “controlled”  transition  experiments  in  the 
Low  Turbulence  Variable  Geometry  wind  tunnel.  Typical  results  from  our 
simulations  of  the  “natural”  transition  experiments  were  presented  in  this 
report.  To  confirm  that  simulations  and  experiments  were  based  on  the 
same  “baseflow,”  the  experimental  baseflow  profiles  were  compared  with  our 
Navier-Stokes  results.  The  agreement  was  reasonably  good,  except  for  the 
last  downstream  location  which  was  outside  the  “quiet”  region.  The  down¬ 
stream  development  and  the  spatial  growth  rates  of  the  disturbances  obtained 
from  the  Navier-Stokes  computations  were  compared  with  the  experimental 
measurements.  The  agreement  was  remarkable,  considering  that  the  data 
were  from  natural  transition  experiments  and  considering  the  difficulties  and 
uncertainties  when  trying  to  carry  out  supersonic  stability  and  transition 
experiments. 

From  these  and  other  comparisons  of  our  Navier-Stokes  results  with  ex¬ 
perimental  data,  we  are  reasonably  sure  that  more  challenging  transition 
experiments  can  be  performed  in  the  experimental  facility  such  as,  for  ex¬ 
ample,  investigations  of  the  nonlinear  stages  of  transition  and  the  onset  of 
the  breakdown  to  turbulence.  Towards  this  end,  G.  Brown  has  been  devel¬ 
oping  a  technique  to  allow  for  “controlled”  transition  experiments  using  a 
loudspeaker  upstream  of  the  contraction  nozzle  to  generate  controlled  (for 
example  single  frequency)  disturbances  (Brown  &  Fan,  2003).  The  experi¬ 
mental  data  obtained  so  far  clearly  indicate  that  “controlled”  experiments 
can  indeed  be  carried  out  in  this  facility.  In  addition,  because  of  the  forc¬ 
ing  technique  employed  in  these  experiments,  the  receptivity  to  sound  can 
also  be  investigated.  In  order  to  perform  such  receptivity  investigations,  we 
have  already  modified  our  code.  Towards  this  end,  we  first  had  to  implement 
free-stream  boundary  conditions  to  allow  for  variable  free-stream  pressure. 
To  validate  the  modified  code,  transition  simulations  with  adverse  pressure- 
gradients  for  a  flat  plate  were  performed  at  low  Mach  number  ( Ma  =  0.1) 
and  the  results  were  compared  to  simulations  from  an  incompressible  code. 
Towards  the  understanding  of  the  nonlinear  mechanisms,  we  simulated  nu¬ 
merous  nonlinear  resonance  and  breakdown  scenarios  for  a  Mach  3  flat-plate 
boundary  layer  with  and  without  adverse  pressure  gradients.  The  objective 
of  these  simulations  was  to  identify  scenarios  that  may  be  relevant  for  labora¬ 
tory  and/or  realistic  flight  conditions.  With  such  simulations  we  can  screen 
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the  extremely  large  parameter  space  for  possible  nonlinear  interaction  and 
breakdown  scenarios  so  that  the  most  viable  and  relevant  ones  can  be  investi¬ 
gated  rigorously  in  experiments  and  detailed  simulations.  This  will  result  in 
considerable  savings  of  time  and  experimental  costs  of  future  investigations. 


APPENDIX  A:  Code  Parallelization 


For  the  present  investigation,  most  simulations  have  been  carried  out  on  ma¬ 
chines  of  the  SGI  Origin  3k  series  with  MIPS  400-700  MHz  Rl2k  or  Rl6k 
processors,  IP35  or  IP27  architecture  and  8  MB  cache  and  HP  Alpha  ma¬ 
chines.  The  efficiency  of  the  parallelization  algorithm  for  these  shared  mem¬ 
ory  machines  has  been  tested  and  it  turned  out  to  be  high  enough  to  justify 
the  use  of  up  to  400  processors  —  the  maximum  number  of  CPUs  possible 
for  the  largest  problem  size  investigated.  For  some  calculations  even  a  su- 
perlinear  Speedup  could  be  obtained,  i.e.,  running  a  simulation  with  n-times 
the  processors  is  more  than  just  n-times  faster.  The  reason  for  this  lies  in 
the  fact  that  for  smaller  computations  the  memory  requirements  can  exceed 
the  cache  of  one  processor  but  do  not  exceed  the  combined  cache  available 
to  the  n  processors.  Thus,  if  feasible,  the  number  of  processors  was  chosen 
to  be  sufficiently  large  for  exploiting  cache  effects. 


In  order  to  evaluate  the  performance  of  the  parallel  algorithm,  we  define  the 
speedup 


Sp  = 


t(n) 


(A.1) 


and  the  parallel  efficiency 


(A.2) 


where  t(n)  is  the  elapsed  time  of  the  computation  using  n  processors.  Am¬ 
dahl’s  law  states  that  the  speedup  is  limited  by  the  parallel  fraction  fp  of  the 
program: 
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for  n  — >  oo.  A  lower  bound  of  the  performance  of  any  program  can  be  esti¬ 
mated  in  millions  of  floating-point  operations  per  second  (MFLOPS). 


Program  performance  on  an  SGI  Origin  3800  with  IP35  architecture  is  evalu¬ 
ated  by  computing  200  timesteps  of  a  three-dimensional  turbulent  boundary 
layer  simulation  using  240  x  130  points  and  16  Fourier  modes  (33  physical 
points).  It  represents  a  typical  calculation  which  includes  data  input  and 
output.  The  number  of  processors  are  varied  and  the  problem  size  is  kept 
constant,  thus  increasing  the  relative  overhead  of  the  parallelization  with  in¬ 
creasing  number  of  CPUs.  The  memory  requirements  for  this  case  are  such 
that  more  than  33  processors  are  needed  to  fit  the  memory  completely  into 
cache,  and  40  processors  are  the  maximum  possible  due  to  domain  length 
limitations. 


Number  of 
CPUs 

MFLOPS 
per  process 

Total 

MFLOPS 

CPU 
time  [s] 

Sp 

Bp 

1 

101 

101 

2,500 

1.0 

100% 

2 

106 

212 

1,197 

2.1 

104% 

4 

94 

376 

691 

3.6 

91% 

8 

90 

720 

348 

7.2 

90% 

12 

89 

1068 

239 

10.5 

87% 

16 

80 

1280 

198  ■ 

12.6 

79% 

20 

72 

1440 

177 

14.1 

71% 

24 

77 

1848 

142 

17.6 

73% 

30 

84 

2520 

106 

23.6 

79% 

36 

87 

3132 

86 

29.1 

81% 

40 

85 

3400 

79 

31.7 

79% 

Table  A.l  Program  performance  on  an  SGI  Origin  3800  for  different  numbers 
of  processors  and  a  fixed  problem  size. 


Table  A.l  and  figure  A.l  clearly  show  that  increasing  the  number  of  proces¬ 
sors  incurs  a  performance  penalty  because  of  the  additional  communication 
associated  with  the  parallelization.  The  penalty  is  visible  in  the  decrease  of 
MFLOPS  per  processor  and  parallel  efficiency.  Nevertheless,  overall  program 
performance  (speedup  and  total  MFLOPS)  is  still  increasing  and  computa¬ 
tion  time  is  decreasing.  Starting  with  the  use  of  20  CPUs,  the  communication 
penalty  is  more  than  compensated  by  the  appearance  of  cache  effects.  With 
increasing  number  of  processors  and  hence  increasing  total  cache  size,  more 
and  more  of  the  data  can  be  accessed  faster  and  parallel  efficiency  increases 
considerably.  At  36  processors  all  data  fit  in  cache.  A  further  increase  in 
CPUs  does  not  enhance  any  cache  effects  and,  consequently,  can  not  compen¬ 
sate  any  concomitant  increase  in  communication  penalty.  Thus  the  efficiency 
decreases  again.  However,  the  theoretical  lines  in  figure  A.l  show  that  fitting 
the  data  in  cache  is  still  equivalent  to  computing  with  an  algorithm  of  higher 
parallel  fraction  —  shifting  the  maximum  speedup  to  higher  values. 

In  addition  to  exploiting  cache  effects,  program  performance  can  be  enhanced 
by  optimizing  the  scheduling  environment  of  the  multi- processor  computer. 
In  figure  A.l,  results  using  the  system  default  environment  (squares)  are 
compared  with  a  user  defined  environment  (solid  circles).  For  more  than  four 
CPUs,  the  effect  of  an  optimized  scheduler  is  again  equivalent  to  increasing 
the  parallel  fraction  of  the  algorithm  and  hence  highly  effective  in  the  case  of 
a  larger  number  of  processors.  Clearly,  any  time  invested  in  optimizing  the 
scheduler  environment  is  time  well  spent. 
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Figure  A.l  Speedup  (left)  and  efficiency  (right)  of  parallization  on  an  SGI 
Origin  3800 ;  circles  (•)  and  squares  (□)  from  computation  with  and  without 
optimized  multiprocessing  environmental  variables,  respectively;  theoretical 

lines  from  Amdahl’s  law  with  parallel  fractions  - —  100%, - 99.3%,  — 

99%,  and  -  •  -  90%. 


APPENDIX  B:  Grid  Convergence  Studies 


B.l  Linear  Stability  Investigations 


Figure  B.l  Influence  of  different  domain  heights  on  amplification  rate  (32 
points  per  Ax,  10  points  per  <5). 


B.1.1  Domain  Height 

Computations  with  3,  6,  9  and  12  6 outflow  are  carried  out  to  assure  that  the 
free-stream  boundary  is  not  influencing  the  results.  All  other  parameters 
stay  the  same  while  the  number  of  points  in  the  y-direction  are  increased. 
Figure  B.l  shows  that  three  boundary  layer  thicknesses  are  too  low  and, 
therefore,  the  amplification  rates  differ,  especially  where  the  wave  is  highly 
amplified,  from  the  other  three  cases.  The  cases  with  six  and  nine  boundary 
layer  thicknesses  in  the  domain  are  identical  except  the  fact  that  the  expan¬ 
sion  fan  originating  from  the  disturbance  slot  hits  the  free-stream  boundary 
later  and  causes  a  slight  oscillation  of  the  amplification  rate  farther  down¬ 
stream.  For  the  three  boundary  layer  thickness  case,  this  modulation  is 
farther  upstream  and  therefore  not  visible  in  figure  B.l. 

The  domain  height  can  also  have  an  influence  on  the  amplitude  and  phase 
distribution.  Figure  B.2  shows  slight  deviations  of  the  free-stream  (77  >  10) 
only  for  a  domain  height  of  3  <W/w 

From  these  results,  a  domain  height  of  6  5^^  should  be  reasonable  for 
further  computations,  but  the  oscillations  in  the  amplification  rate  are  un¬ 
satisfactory.  The  computation  with  12  <W«0U,  is  high  enough  so  that  the 


Figure  B.2  Influence  of  different  domain  heights  on  the  amplitude  distribu¬ 
tions  (left  column)  and  phase  (right  column)  of  u',  ft  and  T  at  Rx  =  500  (32 
points  per  \x,  10  points  per  J). 


expansion  fan  does  not  hit  the  free-stream  and  no  wiggles  are  present.  It 
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Figure  B.3  Influence  of  the  domain  length  (32  points  per  A*,  10  points  per 
5 ,  stretched  9  5  high). 

shows  a  perfect  alignment  with  the  9  5outflow  case  upstream  of  the  oscilla¬ 
tions.  Because  of  the  computational  cost  an  equidistant  grid  is  unreasonable. 
Therefore,  stretched  grid  studies  are  also  performed  (see  appendix  B.1.7). 

B.1.2  Domain  Length 

To  ensure  a  domain  without  any  interferences  of  the  inflow  and  outflow 
boundaries,  the  overall  domain  length  is  under  investigation  in  this  section. 
In  figure  B.3  all  results  are  perfectly  aligned  demonstrating  that  the  upstream 
influence  of  the  outflow  is  negligible.  Figure  B.3  also  shows  that  three  wave¬ 
lengths  after  the  disturbance  slot  (x2)  this  wave  has  not  fully  developed  before 
propagating  out  of  the  domain.  The  five,  seven  and  nine  wavelength  cases 
all  coincide  leading  to  the  conclusion  that  the  domain  length  has  to  be  at 
least  5  A*  long.  Nonetheless,  the  overall  domain  length  also  depends  on  the 
region  of  interest. 

B.1.3  Inflow  Location 

If  the  inflow  boundary  is  too  close  to  the  disturbance  slot,  a  modulation  of 
the  amplification  rate  can  be  experienced.  To  exclude  this  possible  influence, 
four  different  locations  of  the  inflow  are  investigated,  where  the  disturbance 
slot  stays  at  a  fixed  physical  position.  Figure  B.4  shows  that  there  is  no 
visible  difference  between  the  one  wavelength-case  and  the  half  wavelength 
case.  For  a  spacing  of  a  quarter  wavelength  between  the  disturbance  slot 
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Figure  B.4  Influence  of  the  inflow  boundary  location  (32  points  per  A*,  10 
points  per  <5,  9  <5  high). 


and  the  inflow  boundary,  the  amplification  rate  is  first  larger  than  for  the 
other  cases,  but,  from  Rx  =  500  on,  all  amplification  rates  agree  well.  For 
a  distance  of  Ax/8,  the  amplification  rate  curves  do  not  differ  from  the  case 
with  the  largest  largest  spacing. 

B.1.4  Disturbance  Slot  Location 

The  influence  of  the  inflow  on  the  amplification  rate  is  investigated  by  look¬ 
ing  at  different  downstream  locations  of  the  disturbance  slot.  The  chosen 
distances  of  Ax/2,  Ax  and  2AX  show  that  the  location  of  the  disturbance  slot 
does  not  significantly  influence  the  amplification  rates  far  downstream  (see 
figure  B.5).  Farther  upstream  (Rx  <  650),  the  differences  in  the  amplifica¬ 
tion  rates  show  that  there  is  still  some  adjustment  and  wave  development  for 
the  2Ax-case.  Due  to  the  outcome  of  the  study  on  different  inflow  locations 
(see  appendix  B.1.3)  the  differences  between  the  two  cases  with  Ax/2  and  Ax 
stam  from  the  ongoing  adjustment  of  the  wave  of  the  Ax-case.  Because  the 
differences  of  all  three  cases  are  within  acceptable  limits  and  a  farther  down¬ 
stream  location  is  more  applicable  for  this  study,  a  disturbance  slot  location 
of  one  wavelength  downstream  of  the  inflow  is  chosen. 

B.1.5  Outflow  Location 

For  the  linear  stability  investigations  in  6.2,  no  relaminarization  zone  is  used 
because  only  small  amplitudes  are  under  investigation  and  consequently  a 
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Figure  B.5  Influence  of  the  disturbance  slot  location  (32  points  per  Xx,  10 
points  per  6,9  6  high). 


Figure  B.6  Upstream  influence  of  the  outflow  boundary  (32  points  per  A*, 
10  points  per  6,6  6  high). 

buffer  zone  would  only  increase  the  number  of  points  without  enlarging  the 
usable  domain.  Therefore,  the  upstream  influence  of  the  outflow  boundary 
condition  has  to  be  examined.  For  a  frequency  of  F  =  8.1  x  10-5,  the  last 
point  of  interest  is  at  about  Rx  =  800.  From  that  point  on,  the  domain 
is  enlarged  downstream  by  1,  2  and  3  Ax  for  this  investigation.  Figure  B.6 
shows  that  two  wavelengths  are  sufficient  to  avoid  any  upstream  influences. 

B.1.6  Wall-Normal  Resolution 

To  investigate  how  the  resolution  in  the  y-direction  is  influencing  the  ampli¬ 
fication  rates,  computations  with  5,  10  and  15  points  per  6injiam  are  carried 
out.  Figure  B.7  shows  the  low  resolution  case  results  in  a  wrong  amplifica- 
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Figure  B.7  Low  resolution  case  (left)  and  higher  resolution  cases  (right)  in 
the  y-direction  (32  points  per  Ax,  65  high). 


tion  rate.  The  two  higher  resolution  cases  (see  also  figure  B.7)  collapse  to 
one  curve  suggesting  that  ten  points  per  5inf[ow  are  sufficient  to  capture  the 
right  amplification  rate. 

B.1.7  Grid  Stretching  in  the  Wall-Normal  Direction 

Not  to  influence  any  physical  properties  of  the  flow,  the  grid  is  only  stretched 
above  two  boundary  layer  thicknesses.  Thus,  the  second  maximum  of  the 
amplitude  distribution  remains  in  the  equidistant  region.  The  stretching 
function  enlarges  the  domain  to  6  and  9  <50ut/foun  keeping  the  stretching  fac¬ 
tor  constant.  Comparing  the  results  to  the  computations  performed  on  an 
equidistant  grid  show  that  the  wiggles  of  the  amplification  rate  disappeared. 
The  reason  for  this  that  the  grid  stretching  introduced  an  additional  damping 
by  an  artificial  viscosity  caused  by  the  stretching  function.  So  the  expansion 
fan  hitting  the  free-stream  boundary  is  an  order  of  magnitude  smaller  and, 
therefore,  small  enough  not  to  influence  the  amplification  rate  any  more. 
Nevertheless,  the  six  and  nine  domain  height  cases  agree  well  with  the  12- 
5 out  flow  case  on  an  equidistant  grid  (see  figure  B.8).  The  lowest  domain  height 
for  a  stretched  grid  still  shows  a  slight  oscillation  revealing  that  the  reflec¬ 
tion  from  the  shock  wave  still  has  a  small  impact  on  the  amplification  rate. 
The  stretched  9  and  equidistant  12 -Soutfiow  cases  are  in  perfect  agreement 
(see  figure  B.9),  hence  a  domain  height  of  9  boundary  layer  thicknesses  on  a 
stretched  grid  is  suited  for  this  investigation.  To  verify  that  the  stretching 
has  no  influence  on  the  amplitude  and  phase  distribution,  a  comparison  with 
the  equidistant  12  <5out//ourCase  is  performed  (see  figure  B.10). 

The  resulting  domain  height  of  nine  boundary  layer  thicknesses  seems  high, 


Figure  B.8  Comparison  of  stretched  with  equidistant  grid 


Figure  B.9  Comparison  of  stretched  grid  to  equidistant  12  S  (32  points  per 
Ax,  10  points  per  S). 
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Figure  B.10  Influence  of  the  grid  stretching  on  the  amplitude  distribution 
(left  column)  and  phase  (right  column)  of  u\  p'  and  V  at  R,  =  500  (32 
points  per  Xx,  10  points  per  S). 


Ill 


Figure  B.ll  Influence  of  resolution  in  x  on  the  amplification  rate  on  a 
stretched  grid  (10  points  per  5,9  5  high). 

but  as  pointed  out  in  section  3.4.4,  the  flow  is  assumed  to  be  steady  at  the 
free-stream.  For  the  high  frequencies  under  investigation  this  is  a  rough 
approximation,  so  the  domain  height  increases.  The  additional  damping 
caused  by  the  stretched  grid  and  the  fact  that  only  a  third  of  the  points  as 
in  the  equivalent  equidistant  case  are  used,  clearly  shows  that  the  stretched 
grid  is  not  only  the  better  choice  from  a  physical  point  of  view  but  also  from 
a  computational  point  of  view,  reducing  the  computation  time  drastically. 

B.1.8  Streamwise  Resolution 

Computations  are  performed  on  a  grid  which  is  equidistant  over  the  first 
two  boundary  layer  thicknesses  from  the  wall  and  is  then  stretched  to  a 
total  height  of  nine  boundary  layer  thicknesses  (chosen  in  accordance  to  sec¬ 
tion  B.  1.7).  The  grid  in  the  downstream  direction  is  equidistant.  The  number 
of  grid  points  per  TS-wavelength  are  varied  to  determine  the  regular  reso¬ 
lution.  The  16  points  per  TS-wave  case  still  shows  strong  oscillations,  how¬ 
ever,  the  two  higher  resolutions  are  in  very  good  agreement(see  figure  B.ll). 
Therefore,  32  points  per  wavelength  are  chosen  in  the  stability  investigations. 

B.2  Oblique  Breakdown  Simulations 

For  the  validation  cases  and  the  linear  stability  investigations,  all  calculations 
were  performed  on  an  equidistant  grid.  Since  the  resolution  does  not  need 


to  be  high  to  simulate  3D  waves  in  the  linear  stage,  an  equidistant  grid 
is  not  a  disadvantage.  For  oblique  breakdown  this  is  not  true  anymore. 
When  the  3D  waves  start  to  interact  with  each  other,  smaller  structures  are 
being  created  and  these  structures  need  to  be  resolved.  For  that  reason, 
a  high  resolution  is  required  at  positions  where  smaller  structures  emerge 
and  not  further  upstream  where  the  3D  waves  are  in  the  linear  stage.  The 
DNS  code  used  for  all  presented  simulations  has  the  capability  to  use  a 
stretched  grid  in  the  x-  and  in  the  y-direction.  This  feature  was  applied 
for  all  oblique  breakdown  simulations.  Figure  B.12  shows  the  computational 
grid  used  for  the  simulation  of  case  OBLNP2.  The  grids  for  all  other  oblique 
breakdown  simulations  (cases  OBLNP3,  OBLNP4  and  OBLPG)  are  similar. 
The  parameters  which  are  necessary  to  create  such  a  grid  can  be  found  in 


Figure  B.12  Computational  grid  for  the  simulation  OBLNP2.  Every  fifth 
point  is  plotted  for  clarity. 


Figure  B.13  Parameters  for  the  generation  of  the  computational  grid. 


figure  B.13.  For  the  x-direction,  xsi  is  the  location  where  the  stretching 
starts.  The  equidistant  grid  upstream  of  xsi  has  the  stepsize  Ax.  At  position 
xs2  the  stretching  ends  and  downstream  of  xs2  the  grid  is  equidistant  again 
with  the  stepsize  Ax.  The  stepsize  Ax  is  calculated  from  the  stepsize  Ax 


and  the  specified  ratio  Ax /Ax.  For  the  oblique  breakdown  calculations  this 
ratio  is  smaller  than  1,  so  that  the  grid  becomes  finer  with  increasing  x.  For 
the  ?/- direction,  ys\  is  the  location  where  the  stretching  starts.  The  stretching 
ends  at  the  upper  boundary,  that  means  at  x/m  with  the  stepsize  Ay- Ay /Ay. 
Simulation  OBLNP2  is  a  test  case  which  was  conducted  to  find  out  if  the  grid 
stretching  has  an  influence  on  the  amplification  rate  a,  of  a  3D  disturbance  in 
the  iinear  stage.  Case  OBLNP1  is  a  calculation  performed  on  an  equidistant 
grid  with  a  fine  resolution.  Both  calculations  have  the  same  flow  and  forcing 
properties  as  the  simulations  for  the  supersonic  validation.  A  more  detailed 
summary  of  the  input  data  for  the  DNS  can  be  found  in  appendix  C.  The 
amplitude  distribution  over  x  of  the  u- velocity  of  case  OBLNPl  serves  as 
comparison  for  case  OBLNP2.  Since  the  grid  of  case  OBLNPl  is  finer  than  for 
the  supersonic  validation  cases  presented  in  §  4,  the  amplitude  distribution  of 
case  OBLNPl  describes  the  correct  physical  behavior.  Figure  B.14  contains 
the  results  of  both  calculations.  The  slope  of  the  curves  are  the  same  which 
means  that  both  simulations  produce  the  same  amplification  rate  a,. 


Figure  B.14  Amplitude  distribution  versus  x  for  cases  OBLNPl  and 
OBLNP2. 


APPENDIX  C:  Computational  Details 


The  parameters  that  are  used  for  the  computational  simulations  are  listed  in 
the  following  tables.  All  baseflow  calculations  are  performed  with  adiabatic 
walls  (unless  stated  otherwise),  whereas  the  (unsteady)  forced  flow  calcula¬ 
tions  are  carried  out  with  isothermal  walls,  -i.e.,  the  baseflow  temperature 
distribution. 


Thumm  A1 

Eissler  C 

Physical  parameters: 

Ma 

1.6 

4.8 

Re 

100,000 

100,000 

Pr 

0.71 

0.71 

K 

1.4 

1.4 

^ wall 

adiabatic 

650  K 

Twall 

0 

0 

To o 

300  K 

220  K 

Poo 

101.3  kPa 

2.5501  kPa 

Computational  parameters: 

nx 

300 

800 

ny 

150 

121 

At 

0.001256 

0.0062831853 

Aa: 

0.03704 

0.0238 

A Vwall 

0.00125 

0.007 

AVfs 

0.1875 

0.007 

K 

0.4133 

0.5711987 

Xq 

0.225 

44.1 

V stretch 

0.125 

N/A 

Vm 

3.0 

0.847 

X>ramp 

N/A 

51.5 

Forcing  parameters: 

0 

5.0025 

10 

Ai,o 

0 

5  x  10~2 

^1,1 

5  x  10"5 

1  x  10~3 

Xi 

1.44732 

44.9568 

X2 

1.8332 

45.2424 

Table  C.l  Parameters  for  the  supersonic  validation  cases  of  Thumm  (1991) 
and  Eissler  (1995). 


ICVALDC 

ICVALFR 

CVALNP 

CVALPG 

Physical  parameters: 

Re 

Ma 

T* 

Pr 

K 

pH 

105 

0.25 

300AT 

0.71 

1.4 

0.0 

105 

0.1 

224A 

0.71 

1.4 

-0.18 

1578102 

3.0 

103.6  AT 

0.71 

1.4 

0.0 

1578102 

3.0 

103.6A: 

0.71 

1.4 

-0.04...  -  0.08 

Domain  size: 

Xq 

Vm 

K 

0.54 

4.124 

«0.1 

N/A 

1.585 

5.844 

0.22 

2.0520  x  10"1 

0.0396 

2.2716 

0.0199 

4.4239  x  10“2 

0.5116  (0.0396) 
2.1036  (3.1920) 
0.0199 

3.7920  x  10“2 

Buffer  domains  and  decay  condition: 

Jin 
^  start 

rpin 

^ end 
~out 
'''start 

rpOUt 

■ end 

a 

C 

N/A 

N/A 

N/A 

N/A 

37.6 

N/A 

1.673 

2.060 

2.800 

4.200 

29.8 

N/A 

N/A 

N/A 

1.8796 

2.1756 

70.0 

0.7 

N/A 

N/A 

1.7036 

2.0236 

74.0 

0.7 

Grid  size  and  resolution: 

nx 

ny 

K 

At 

Ax 

Ay 

513 

88 

0 

1  x  10~4 

7  x  10‘3 

5  x  10"4 

506 

292 

2 

4.8481  x  10-5 
8.4340  x  10"3 
5.7300  x  10"4 

280 

200 

1 

3.3179  x  10"5 
8.0000  x  10~3 
1.0000  x  10"4 

200  (400) 

200 

1 

3.3179  x  10~5 
8.0000  x  10"3 
1.0000  x  10"4 

Forcing  parameters: 

P 

^-1,0 

Xi 

%2 

14.0 

1  x  10"4 
N/A 

0.82 

0.89 

10.8 

1  x  10“6 

2  x  10~6 
2.082 

2.293 

47.343 

N/A 

1  x  10“6 
0.2236 
0.2716 

47.343 

N/A 

1  x  10"12 
0.6956 

0.7436 

Table  C.2  Parameters  for  the  pressure-gradient  validation  cases;  baseflow 
computations  are  carried  out  in  2D,  values  which  differ  otherwise  are  given 
in  brackets;  pressure-gradient  parameters  for  case  CVALPG  are:  x$  —  0.6, 
Xi  =  2.0,  xc,  —  2.6,  Api  =  0.8,  A p2  =  0.6  and  transdx  —  0.1. 


F  =  1.4  •  10~5 

|  .F  =  5.0  •  10~5 

F  -  8.1  •  10"5 

Physical  parameters: 

Ma 

3.0 

3.0 

3.0 

Re 

1,578,102 

1,578,102 

1,578,102 

Pr 

0.71 

0.71 

0.71 

K 

1.4 

1.4 

1.4 

Twall 

adiabatic 

adiabatic 

adiabatic 

V 

wall 

0 

0 

0 

Too 

103.6  K 

103.6  K 

103.6  K 

Poo 

0.7508  kPa 

0.7508  kPa 

0.7508  kPa 

Computational  parameters: 

'ft'X 

400 

280 

434 

Uy 

236 

328 

182 

At 

5.636138  •  10"5 

1.61085  •  10~5 

4.973063  •  10~5 

Ax 

4.20175  •  10~3 

1.6228 -lO"3 

1.07456  •  10"3 

A Vwall 

9.72584  •  10-5 

3.51  •  lO”5 

6.3776 -lO"5 

&Vjs 

9.72584  •  10~5 

2.2113  •  10"3 

1.721952  •  10“3 

A* 

0.0970526 

0.0299817 

0.0198527 

x0 

0.017836 

1.591725  •  lO'3 

6.3669 -lO"3 

y  stretch 

N/A 

9.7578  •  10"3 

8.3706 -lO"3 

Vm 

0.022953 

0.047 

0.0392 

Forcing  parameters: 

/? 

22.29607 

•  78.0107 

126.34437 

AU,i 

5 • 10"4 

5 • 10~5 

5  •  10“4 

Xl 

0.08086225 

0.027556525 

0.04075282 

X2 

0.101871 

0.040519725 

0.05794578 

Table  C.3  Parameters  for  the  simulations  for  comparison  with  the  experi¬ 
ments  of  Graziosi  (1999)  and  Graziosi  &  Brown  (2002). 
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II J  1  H'BM'MI 

Physical  parameters: 

Ma 

3.0 

3.0 

3.0 

1  .  3.0 

Re 

1,578,102 

1,578,102 

1,578,102 

1,578,102 

Pr 

0.71 

0.71 

0.71 

0.71 

K 

1.4 

1.4 

1.4 

1.4 

TWaa 

adiabatic 

adiabatic 

adiabatic 

adiabatic 

mr 

wall 

0 

0 

0 

0 

103.6  K 

103.6  K 

103.6  K 

103.6  K 

Poo 

0.7508  kPa 

0.7508  kPa 

0.7508  kPa 

0.7508  kPa 

Computational  parameters: 

nx 

600 

570 

384 

384 

Uy 

174 

134 

180 

224 

At 

8.05426  •  10~5 

4.973063  •  10~5 

2.651  •  10~s 

1.9899  •  10"s 

Ax 

1.6228  •  10~3 

1.07456  •  10-3 

7.3125  •  10~4 

4.63  •  10"4 

^ Vwall 

1.3299  -10-4 

1.40307-  10~4 

6.63  •  10~5 

3.536  •  10~5 

7.9794  •  10“3 

2.806144  •  10“3 

4.641  •  10~3 

6.3648  •  10~3 

A  , 

0.0299817 

0.0198527 

0.01351 

8.55384  •  10"3 

Xo 

0.3105 

0.0254676 

1.591725  •  10"3 

1.591725  •  10"3 

y stretch 

0.0170227 

0.017858 

8.619  •  lO"3 

6.15264  •  lO"3 

Vm 

0.143 

0.0631 

0.087 

0.111 

Forcing  parameters: 

P 

78.0107 

126.34437 

237 

315.758 

Au 

5 • 10-5 

5 • 10"4 

5  •  lO"5 

5 • 10~5 

*1 

0.3624296 

0.05985352 

0.024991725 

8.9997  •  10~5 

0.3883944 

0.07704648 

0.036691725 

0.01641325 

Table  C.4  Parameters  for  the  simulations  for  investigation  of  non-parallel 
effects. 


OBLNP1 

OBLNP2 

OBLNP3 

OBLNP4 

OBLPG 

Physical  parameters: 

Re 

1578102 

|  1578102 

1578102 

1578102 

1578102 

Ma 

3.0 

3.0 

3.0 

3.0 

3.0 

rp* 

1  OO 

103.6iir 

103.6iT 

103.6# 

103.6# 

103.6# 

Pr 

0.71 

0.71 

0.71 

0.71' 

0.71 

K 

1.4 

1.4 

1.4 

1.4 

1.4 

Ph 

0.0 

0.0 

0.0 

0.0 

«  -0.12 

Domain  size: 

Xo 

0.35644 

0.35644 

0.35644 

0.95244 

0.95244  (0.35644) 

xN 

1.31244 

1.31211 

1.3095 

2.01177 

2.01177  (2.48277) 

Vm 

0.0369 

0.03687 

0.0366 

0.037 

0.037 

K 

4.10  x  10~2 

4.10  x  IQ"2 

4.10  x  lO"2 

4.10  x  lO"2 

3.73  x  lO"2 

Buffer  domains  and  decay  condition: 

™out 
^  start 

N/A 

N/A 

1.20644 

1.86277 

1.86277 

rpOUt 

^ end 

N/A 

N/A 

1.30644 

2.00277 

2.00277 

a 

70.0 

70.0 

70.0 

70.0 

78.5 

c 

0.7 

0.7 

.0.7 

0.7 

0.7 

Grid  size,  resolution  and  stretchinj 

7- 

nx 

240 

430 

585 

550 

550  (1170) 

ny 

370 

80 

130 

140 

140 

K 

1 

1 

.  54 

15 

15 

At 

1.33  x  10~4 

1.33  x  lO'4 

1.33  x  10-4 

1.33  x  10"4 

1.33  x  10~4 

Ax 

4  x  10"3 

4  x  10-3 

4  x  10~3 

4  x  10-3 

4  x  10“3 

Ay 

1  x  10"4 

1  x  10"4 

1 

O 

i — 1 

X 

iH 

1  x  10"4 

1  x  10"4 

XS1 

N/A 

0.63244 

0.63244 

1.25644 

1.25644 

XS2 

N/A 

1.18044 

1.18044 

1.74044 

1.74044 

Ax 

Ax 

N/A 

0.25 

0.125 

0.25 

0.25 

ys  1 

N/A 

0.0 

0.005 

0.0 

0.0 

*SL 

Ay 

N/A 

12.0 

10.0 

6.0 

6.0 

Forcing  parameters: 

P 

47.343 

47.343 

47.343 

47.343 

47.343 

^1,1 

10~5 

O 

1 

01 

0.003 

0.0054 

0.0014 

Xl 

0.54444 

0.54444 

0.54444 

1.12444 

1.12444 

X2 

0.62444 

0.62444 

0.62444 

1.19644 

1.19644 

Table  C.5  Parameters  for  the  oblique  breakdown  simulations;  baseflow  com¬ 
putations  are  carried  out  in  2D,  values  which  differ  otherwise  are  given  in 
brackets;  pressure-gradient  parameters  for  case  OBLPG  are:  x3  =  1.15644, 
£4  =  1.97244,  £5  =  2.27644,  Api  =  0.53,  A p2  =  0.53  and  transdx  =  0.1. 
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